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ABSTRACT 


In  solving  the  Boltzmann  transport  equation,  most  discrete  ordinates  codes 
calculate  the  source  term  by  first  approximating  the  scattering  cross  section 
using  a  Legendre  polynomial  expansion.  Such  expansions  are  insufficient  when 
scattering  is  anisotropic  and  the  Legendre  expansion  is  truncated  prematurely. 
This  can  lead  to  nonphysical  negative  cross  sections,  negative  source  terms  and 
negative  angular  fluxes.  While  negative  sources  are  problematic  for  standard 
discrete  ordinates  methods  leading  to  poor  convergence  or  convergence  to 
incorrect  results,  they  are  of  particular  concern  to  exponential  methods,  causing 
such  calculations  to  fail. 

We’ve  developed  and  tested  a  new  technique  to  solve  this  problem  called  the 
Monte  Carlo  Facet  Method.  This  method  is  an  extension  of  standard  Monte 
Carlo  techniques.  It  guarantees  non-negative  cross  sections  at  all  directional 
ordinates.  It  also  ensures  within  group  and  next  group  scatter. 

This  dissertation  outlines  previous  attempts  to  handle  anisotropic  scattering 
to  achieve  non-negative  sources.  It  develops  the  theory  of  the  Monte  Carlo  facet 
method  and  its  first  angular  moment  conservation.  Results  are  presented 
examining  the  scattering  matrices  for  various  materials,  and  finally 
demonstrating  that  these  scattering  matrices  perform  exceptionally  well  in  a 
multi-group,  anisotropic,  unstructured  mesh  discrete  ordinates  transport  code. 
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NEUTRON  TRANSPORT 
WITH  HIGHLY  ANISOTROPIC  SCATTERING 


1.  Introduction 

Open  any  standard  transport  textbook  from  the  50's  or  60's,  and  in  the 
chapter  on  discrete  ordinates,  one  of  the  first  things  you'll  need  to  do  is 
approximate  the  scattering  source  with  a  Legendre  polynomial.  These  codes 
were  run  on  machines  where  memory  was  at  a  premium.  The  problem, 
however,  is  that  such  expansions,  among  other  things,  lead  directly  to  the 
calculation  of  negative  scattering  sources  and,  in  turn,  negative  fluxes.  This 
happens  when  the  group-to-group  scattering  cross  section  is  anisotropic. 
Polynomials  can't  easily  approximate  such  behavior  (Figure  1). 


Figure  1.  Within  group  (13.5-14.9  MeV)  elastic  scattering  cross  section  of 

238U  (ref.  5). 
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Today,  we  have  desktop  machines  with  10  gigabyte  hard  drives,  lOO's  of 
megabytes  of  RAM  and  processors  that  outperform  the  mainframes  of 
yesteryear.  Yet,  as  it  was  40  years  ago,  open  up  a  new  transport  textbook, 
and  in  the  chapter  on  discrete  ordinates,  you'll  still  need  to  approximate  the 
scattering  source  with  a  Legendre  polynomial. 

The  point  is  that  even  with  the  advent  of  new  technology,  the 
computational  paradigms  of  the  past  remain.  This  research  breaks  those 
paradigms  and  develops  a  wholly  new  approach  to  create  scattering  cross 
sections  for  discrete  ordinates  codes.  In  so  doing,  we  demonstrate  the  abihty 
to  perform  accurate  deep  penetration,  multi-group,  anisotropic  scattering 
transport  problems  without  unphysical  artifacts.  The  solution  of  such 
problems  is  essential  to  understand  and  model  the  penetration  of  radiation 
through  shields  on  everything  from  nuclear  reactors  to  bomb  shelters  to 
waste  transportation  containers. 

The  dissertation  is  divided  into  a  number  of  chapters.  The  first  few  look 
at  the  discrete  ordinates  method  and  solutions  to  the  scattering  source.  The 
discrete  ordinates  method  is  introduced  and  the  scattering  source  is  defined. 
Various  approaches  used  in  calculating  this  source  for  discrete  ordinates 
transport  calculations  are  then  compared.  In  so  doing,  we  show  how 
spherical  harmonic  (SH)  methods  using  Legendre  polynomials  lead  to 
negative  sources,  and  how  current  attempts  to  solve  such  negativity  are 
inadequate. 
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The  next  few  chapters  develop  the  theory  of  our  proposed  solution,  the 
Monte  Carlo  (MC)  facet  method.  This  method  takes  advantage  of  today’s 
computing  power,  offers  a  direct  link  to  the  physics  of  the  scattering,  and 
guarantees  a  non-negative  scattering  source.  The  theoretical  development  of 
the  method  is  presented  as  well  as  the  general  algorithm  the  code  follows  to 
create  the  non-negative  scattering  matrices.  Also  presented  is  a  method  to 
conserve  the  first  angular  moment  of  the  scatter  (important  in  diffusion  like 
problems)  for  these  new  scattering  matrices. 

The  next  chapter  discusses  the  codes  and  cross  section  libraries  used  for 
this  research.  This  chapter  outlines  the  benchmarking  used  and  explains  the 
basic  operation  of  the  code.  It  also  discusses  the  capabilities  and  limitations 
of  the  codes  developed  as  part  of  this  research. 

The  final  two  chapters  present  the  results.  These  are  broken  into  two 
main  categories: 

1.  Results  examining  the  new  scattering  matrices. 

2.  Results  using  those  matrices  in  a  transport  code. 

The  first  of  these  chapters  demonstrates  the  non-negativity  of  the  new 
approach  and  compares  and  contrasts  this  method  with  other  attempts  to 
solve  the  same  problem.  The  second  chapter  demonstrates  that  we’re  able  to 
perform  real  transport  problems  with  real  materials.  It  also  compares 
transport  results  using  the  facet  method  to  those  using  the  standard 
spherical  harmonic  method.  What  is  found  is  that  the  MC  facet  method 
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provides  superior  performance  while  guaranteeing  non-negative  cross-section 
data. 
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2.  The  Discrete  Ordinates  Method 


This  chapter  outlines  the  discretization  of  the  transport  equation  in 
energy,  angle  and  space.  It  begins  with  the  time  independent,  non¬ 
multiplying,  transport  equation^ 


Jq- V  +  cr(r,£^)j  tf/(r,Q,E) 

=  (r,Q,E)  +  J dE’^ dO.'Gs (r,E'  ->  E,Q'  ■  Q)  y/(jr ,€l’ ,E'), 


(1) 


where  y/  is  the  angular  flux,  <Tis  the  total  macroscopic  neutron  scattering 
cross  section,  and  qin  is  the  time  independent  intrinsic  source.  The  source, 
qij^(r,Q,E)dV dQ.dE ,  is  the  rate  of  source  particles  emitted  in  dV about  r  , 

traveling  in  a  cone  of  directions  dQ  about  Q  with  energies  between  E  and 
E  +  6.E,  and  as  is  the  differential  scattering  cross  section.  In  general. 


CTgfr ,E'  — >  E ,Q'  •  O)  —  a giv ,E’^fg{r ,E'  — ^  E,Q’  — >  O) 

where  as(r,E')  is  the  total  macroscopic  scattering  cross  section  (m-i)  and 
fs(r,E'  ->  E,Q'  Q)dQdE  is  the  conditional  probability  that,  given  that  an 
incident  neutron  of  direction  Q'  and  energy  E'  is  involved  in  a  scattering 
collision,  a  scattered  neutron  will  emerge  from  the  collision  in  the  direction 
interval  dQ  about  Q  with  energy  between  E  and  E+dE.^  The  probability 

density  function  is  normalized  so  that  JJdQdE  fs{r,E'  ->  E,Q'  ->  Q)  equals 

the  expectation  value  of  the  number  of  neutrons  that  emerge  from  such  a 
scattering  collision. 
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At  this  point,  the  scattering  source  is  defined  as 


q^(f,Q,E)  =  J  dE'j  dQ'a,ir,E'  E,Cl'  ■  Q)  i//(r,Q',E') ,  (3) 

where  the  integration  over  on  the  unit  sphere,  %  is  normalized  so  that 

1.  •§  2^  T 

jdn=lfjf  =  i,  (4) 

-1^  0 

and  may  be  approximated  by 

ldQf(Q.)«Y.^JiQj,  (5) 

for  some  quadrature  set,  |(u;„ ,Q„)|  7i  =  1,..,A^| ,  where  the  quadrature  weights 
are  normalized  such  that 

(6) 


Energy  Discretization 

The  BTE,  equation  (1),  is  integrated  over  energy  where  the  energy 
integral  is  partitioned  into  groups.  Namely, 

”  G  G 

|d£:=£  JdE=2jdE.  (7) 

0  ^=1  Eg  g=l  g 
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where  G  is  the  number  of  energy  groups  and  g  is  the  group  number.  In 
general,  Eo>Ei>E2...>Eg=0.  The  first  assumption  in  discrete  ordinates  is 
energy  separability  of  the  angular  flux,  so  that 

y/{f,E,h)^Wg{E)ii/g{r,h),  (8) 


where 

\^dEWg{E)  =  l, 


(9) 


and  Wg  is  constructed  from  an  assumed  energy  distribution,  WiE^  as 


WAE)  = 


W{E) 


\dE  W(E) 


(10) 


A  typical  choice  is  W{E)  -  E  ^  often  called  the  slowing  down  spectrum. 
Here,  the  angular  flux  i//g(r,Q)  is  a  group,  angular  flux,  integrated  over  the 


given  group,  and  remaining  a  distribution  in  Q . 

With  this  assumption,  the  group  form  of  equation  (1)  is^ 

[Q-V  +  o-^(r)]  Wg(rA)  =  qg(rA)  +  qi(r,Q),  (H) 


where  the  group  variables  are  defined  as 


qg(r,Q)=  ^  , 

(12) 

g'=i 

crg(r)  =  l^dEc7ir,E)Wg(E), 

(13) 

q^g{rA)  =  \^dE  q^"'(jr,Q,E), 

(14) 
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and  •  Q)  =  \dE\dE'  (j^(x,E'  ->  E,Q'  •  Q)W^,  {E'). 


This  research  seeks  an  effective  approximation  of  the  group  scattering 
source  term  in  equation  (11),  or  using  equations  (12)  through  (15), 


qg(r,Q)=  Yj  ldCl'jdEjdE'cT\r,E' ^E,Ci'  Q)Wg.(E')i//g.(r,h').  (16) 

g'=i  ^ 


Angular  Discretization 


Standard  discrete  ordinates  theory  assumes  that  equation  (11)  holds  for  N 
distinct  angles  (7i=l,2,..N)  where  an  appropriate  angular  quadrature  is 


applied  giving 


•V  +  <T^(r)]  y/g^{r)^q2n{r)  +  ql^{r) 


where 


Qg)iir)  =  qg{r,an) 


^\da'\dE\dE'a^{f,E'  E,^'-^^)Wg\E')ii/g.{r,a’). 


Applying  the  angular  quadrature,  equation  (5), 


G  N 


g'=ln'==l 


SO  that,  using  the  definition  of  cr®,„ , 
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With  this  notation,  ij/g^n^r) ,  is  the  angular  flux  integrated  over  group  g, 
evaluated  at  the  directional  ordinate Q„ .  Similarly,  Qg^ir)  is  the  scattered 
group  source  evaluated  at  the  directional  ordinate  . 

Material  Dependence 

While  the  source  and  angular  flux  may  vary  through  a  single  material  as 
r  changes,  it  is  assumed  here  that  the  total  and  scattering  cross  sections  do 
not.  While  the  cross  sections  differ  among  materials,  there  is  no  variation 

within  any  single  material.  Hence,  cr^(r,E'  -0,^)  can  really  be 

thought  of  as  E^n'  -^/Jor  cr^^^^iE'  ^  E,Qn'  -^n)  •  Unless 

specifically  required,  the  material  dependence  can  be  assumed  for  all  cross 
sections,  and  the  r  or  material  designation  will  be  left  out  of  further 
derivations.  The  scattering  source  is  then  written  as 

G  N  „  „  ~ 

j^EjdE'a^(E'  ->■  E,Qj^.  ■Q^)Wg>{E')i// g'^,(r)  (21) 

g’=ln'=l  ^  ^ 

Most  modern  transport  codes  approximate  this  scattering  source  with  a 
truncated  spherical  harmonics  expansion.  This  leads  to  the  calculation  of 
negative  sotirces  when  the  scattering  is  anisotropic.  Such  an  expansion,  and 


attempts  to  eliminate  negative  cross  sections  and  hence,  negative  sources  are 
discussed  in  the  following  chapter. 
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3.  Calculating  the  Mulit-Group  Scattering  Cross  Section 


The  previous  chapter  developed  the  standard  discrete  ordinates  theory 
leading  to  the  definition  of  the  scattering  source  term.  This  chapter  examines 
the  standard  method  used  to  calculate  the  scattering  source,  namely  the  use 

of  Legendre  polynomials  to  approximate  (jg’g(Qn’  ’^n)  • 

In  the  past,  direct  calculation  of  the  scattering  source  in  equation  (19)  was 
impossible  because  of  the  memory  requirement  to  store  the  scattering  matrix. 
A  typically  used,  Ss,  quadrature  in  three  dimensions  consists  of  80  angular 
directions.  If  equation  (19)  is  used  directly  for  each  ordinate,  then 

)  is  an  80x80  matrix  for  each  group -to-group  transfer  for  each 
material  used  in  the  problem.  For  a  group  structure  consisting  of  30  groups 
the  cross-section  matrix,  cr|r^(Q„.  •^„)  consists  of  =  5,760,000 

elements.  Stored  as  four-byte  floating-point  numbers,  the  file  containing  the 
matrix  would  be  large,  about  20  Mbytes.  A  further  memory  penalty  is  that 
the  angular  flux  values  for  each  group  and  angle,  y/ ,  must  be  stored  for 

the  source  calculations  of  all  lower  groups  (g>g'),  assuming  down  scatter 
only.  Historically,  authors  of  production  codes  have  considered  this  too  heavy 
a  computational  cost.  With  the  increased  memory  and  storage  capabilities  of 
modern  machines,  the  storage  difficulties  become  less  and  less  an  issue.  In 
fact,  we  will  demonstrate  later  that  such  codes  can  be  run  on  even  a  desktop 
machine.  Still,  problems  requiring  a  large  number  of  materials  and  fine 
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discretizations  in  space,  angle  and  energy  may  not  be  able  to  accept  the 
storage  and  memory  burden  of  direct  calculation.  They  require  that  the 
source  be  approximated  with  spherical  harmonics. 

Using  Spherical  Harmonics 

This  method  is  the  dominant  method  among  production  codes  and 
transport  texts.  The  group-to-group  scattering  cross  section  is  approximated 
byi 

L 

1=0 

where  P;  are  the  Legendre  polynomials  and  cr^ig'g  are  the  Legendre 
coefficients  of  the  group-to-group  scattering  cross  section  found  from 

cj\gg  =  \dE'\dE  af(E'  ^  E)W(E') .  (23) 

g'  g 

The  addition  theorem  for  Legendre  polynomials  is 

p,(n„..Q„)=^  Z  Y:Uh,.)Y,„(h„) ,  (24) 

where  the  Y)t,„  are  the  normalized  spherical  harmonics^  Using  this 
substitution,  the  scattering  cross  section  is 

-QJ «  E  j:  Yl  J •  (25) 

l=Q  m=-l 
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Finally  the  scatter  source  is  found  by  substituting  equation  (25)  into  the 
source  of  equation  (19)  resulting  in 

^gA  *  E  Z  ^n’  E  E  (^n  )Ylm  i^n’  )<^lgg'  ¥ g  fi'  (^)  •  (26) 

^'=l7i'=l  l=0rn=-l 

One  crucial  aspect  in  solving  these  equations  is  that  the  scattering 
coefficients  crigg  are  known  calculated  a  priori  using  codes  such  as  AMPX^  or 

NJOY^.  Still,  once  they  are  known,  for  the  small  cost  of  calculating  the 
spherical  harmonics,  the  memory  requirement  is  significantly  reduced. 

In  a  standard  transport  code,  the  source  is  typically  found  by  first 
calculating  the  moments  of  the  total  group  flux 

<l>Tg  (r) «  E  ^n'Yhn  i^n' )  ¥gg  (r) .  (27) 

n'=l 

Once  moments  of  the  group  flux  are  found,  there  is  no  longer  a  need  for  the 
¥gg  reducing  memory  requirements  further.  The  source  is  found  from 

9j,«(^)“E  E  ^Zm(^«)E  ^igg' <t>\ng(r)  ■  (28) 

l=0m--l  g'=l 

Of  course,  the  spherical  harmonic  approximation  is  just  that,  an 
approximation.  It  is  an  attempt  to  model  the  physics  with  a  polynomial.  In 
many  cases  such  an  approximation  may  be  sufficient.  Unfortunately,  the 
method  becomes  woefully  inadequate  in  a  number  of  ways  particularly  for 
anisotropic  scattering. 
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Anisotropic  scattering  becomes  significant  when  considering 

a.  scattering  of  high  energy  particles  from  nuclei  of  all  mass  numbers, 

b.  scattering  of  any  energy  particle  from  light  nuclei,  and 

c.  when  the  group  structure  is  so  fine  that  the  allowed  angular  region  of 
scatter  is  also  small®. 

Such  scattering  can  be  highly  peaked  and,  because  of  the  group  structure 
used,  be  positive  for  only  a  small  portion  of  the  scattering  range  and  zero  over 
the  rest  of  the  angular  domain. 

Because  ->  Q,j)  =  -Q^),  we  define  =  Q„.  ,  the 

cosine  of  the  scattering  angle  in  the  laboratory  frame.  It  is  often  convenient 
to  examine  the  behavior  of  the  cross-section  in  terms  of  •  We  use 

this  form  of  the  scattering  cross  section  to  demonstrate  the  problems 
Legendre  polynomials  have  approximating  the  type  of  anisotropic  behavior 
described  above.  Figure  2  shows  an  exact  group-to-group  scatter,  , 

for  hydrogen.  Also  shown  in  the  figure  are  three  Legendre  approximations  of 
the  cross  section  using  the  Legendre  expansion  of  equation  (22)  for  L  =  1,  3, 
and  7. 
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Figure  2.  Group  l-»5  scatter  for  hydrogen,  ^H,  in  the  lab  frame. 

The  Legendre  expansion  of  such  a  scatter  fails  in  two  ways.  First,  because 
the  polynomial  is  unable  to  model  such  a  step  function,  it  spreads  the 
distribution  of  scatters  to  other  angles,  causing  an  angular  diffusion  of  the 
scatters  that  is  a  computational  artifact.  Second,  the  expansion  may  become 
negative  in  some  angular  regions.  This  may  lead  to  the  calculation  of 
negative  sources,  which  are  non-physical.  Worse,  the  calculation  of  such 
sources  can  cause  certain  spatial  quadrature  techniques  used  in  transport 
codes  to  faiF-'^.  Particularly  vulnerable  are  exponential  methods  that 
guarantee  non-negative  flux  solutions,  but  require  non-negative  sources®.  In 
the  case  when  negative  sources  appear  in  codes  that  don't  care  or  when  such 
sources  are  fixed  by  setting  them  to  zero  significant  convergence  problems 
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may  arise®.  Such  difficulties  have  led  researchers  to  attempt  alternative 
approaches  that  would  eliminate  negative  cross-section  values. 
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4.  Attempts  to  Fix  Negativity 


A  number  of  attempts  have  been  made  to  eliminate  negative  sources. 

Each  of  the  techniques  described  here  use  the  full  scattering  matrices  to 
calculate  the  scattering  source  term  (possible  now  with  today’s  computational 
capabilities).  This  chapter  outlines  three  techniques,  each  with  its 
shortcomings,  and  introduces  the  method  we  propose. 


Exact  Ordinate-to-Ordinate  Scatter 

Odom  was  one  of  the  first  to  try  to  eliminate  the  problem  of  negative 
sources  in  anisotropic  scattering^®.  To  calculate  the  source  from  equation  (20) 
directly,  he  calculated  the  ordinate-to-ordinate  (Q,i.  Q,j)  scattering  cross 
section  from  the  physics,  but  only  for  elastic  and  level  inelastic  scattering. 
With  this  approach,  he  was  able  to  devise  a  number  of  different  cases  to 
calculate  the  scattering  cross  section  dependent  on  the  angle  and  energies 
involved.  In  particular,  the  group-to-group  scattering  cross  section 


\dE\dE'  o-®(r,E'->E,Q„,  •Q„)W(E') 

/T®  /"O  -  O  \  •'g' 


jdE'  W(E') 


(29) 


was  reduced  to  the  calculation  of  a  1-D  integral®. 

The  integration  must  be  performed  for  each  ordinate-to-ordinate 
combination  (dependent  on  the  quadrature  used),  and  for  each  group-to-group 
combination  (dependent  on  the  group  structure  used).  For  some  scattering 
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combinations,  dependent  on  the  energy-angle  combination,  the  integration 
must  be  performed  more  than  once.  If  level  inelastic  scattering  is  to  be 
calculated,  the  integration  must  be  performed  again  for  each  elastic  and  level 
inelastic  transfer.  In  oxygen,  for  example  there  are  38  inelastic  levels  from 
the  evaluated  nuclear  data  files,  ENDF/B-VI.  The  results  are  then  summed 
to  provide  the  total  transfer  cross-section.  Fortunately  the  integration  is  only 
one  dimensional,  but  the  number  of  cases  that  must  be  considered  make  the 
method  computationally  burdensome,  although  no  mention  is  made  of  the 
time  required  to  compute  the  integration  even  from  more  recent  authors^. 
While  elastic  and  level-inelastic  scatter  can  be  treated  with  a  complicated 
case  structure,  other  types  of  scatter  (e.g.,  (n,  2n))  become  too  difficult  to  put  . 
the  1-D  integral  equation  of  Odom  in  a  readily  integrable  form  and  are 
instead  treated  as  isotropic  even  when  they  are  not. 

While  the  technique  requires  a  complicated  number  of  integration  cases  to 
account  for  variations  in  group  structure,  the  process  guarantees  a  non¬ 
negative  scattering  cross-section.  In  addition  to  the  complexity  of  the 
integration,  the  method  suffers  from  one  major  drawback.  As  Brockman^ 
notes,  for  elements  with  a  highly  anisotropic  angular  distribution  or  for  light 
elements  along  with  a  fine  energy  group  structure,  the  group-to-group 
transfer  cross  section  is  confined  to  a  small  angular  range  (Figure  3).  If  the 
angular  quadrature  used  is  too  sparse  (lack  of  angular  support),  particles 
traveling  in  one  direction  may  never  scatter  into  another. 
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Figure  3.  Group -to-group  Scattering  Geometry^. 


Angular  Support 

Lack  of  angular  support  is  shown  by  examining  the  lab  frame  scattering 
cosine 


ML=S(E,E')^j 


(A- 


(30) 


For  elastic  scattering,  Q  =  0,  and  equation  (30)  can  be  re-written  as 


^JE'-AE  E 

(^+1)1—;:;^ - (^-1) 


E' 


E'-AE 


(31) 


where  AE  is  the  energy  lost  by  the  neutron  in  the  collision.  Expanding 
equation  (31)  in  a  power  series  about  AE  =  0 ,  the  scattering  angle  for  small 
energy  loss  is 

/'L=1-|-^  +  0(aE2).  (32) 


19 


Because,  /ij^  =  cos(0^) » +  small  0i ,  equation  (32)  can  be 


written  in  terms  of  the  scattering  angle 


Ol  a£; 
^  \E' 


(33) 


Using  fine  energy  group  structures  the  within  group  and  next  group 
scattering  angles  will  be  small.  For  within  group  scatter,  the  greatest 
scattering  angle  (smallest  //£,)  occurs  when  a  neutron  scatters  from  the  top  of 
the  group  to  the  bottom  of  the  group  or  h.E  =  where  is  an 

arbitrarily  small  number  greater  than  zero. 

For  ordinate-to-ordinate  methods,  the  6^  ’s  corresponding  to  the  various 
combinations  of  n  and  n'  are  fixed  based  on  the  selected  angular  quadrature. 

To  guarantee  within  group  scatter,  two  directional  ordinates  must  be  close 
enough  so  that 

Sl  (n„.  ^  n„ )  £  ^^(ABg+ofi) .  (34) 


If  is  the  smallest  angle  of  all  the  possible  ordinate-to-ordinate 
scatters,  then  as  the  group  structure  becomes  fine  enough,  AEg  a  2  (where 

again  a2  is  an  arbitrarily  small  number  greater  than  zero),  we  can  find  an 
02  such  that 


0 


min 


(35) 
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and  within  group  scatter  is  not  possible  for  any  of  the  ordinate-to-ordinate 
scatter  combinations. 

The  same  analysis  holds  for  downscattering  to  the  next  lowest  group. 

Here,  the  greatest  scattering  angle  (smallest  fd£)  occurs  when  a  neutron 
scatters  from  the  top  of  the  incident  group  to  the  bottom  of  the  scattered 
group  or  ^E  =  t^Eg  +  A£J^+i  +  ai .  This  scattering  is  not  as  restrictive  on  the 

angular  quadrature  as  within  group  scattering  Because  ->  a2  and 
->  «3  results  in 

<5’min  +«2 +«3)  • 

Still,  it  will  always  be  possible  to  find  a  fine  enough  group  structure  such 
that  downscatter  to  the  next  lower  group  is  not  supported  by  the  angular 
quadrature. 

Does  such  a  failure  in  angular  support  occur  with  typical  angular 
quadratures  and  standard  group  structures?  The  answer  is  yes.  As  is  shown 
in  detail  in  chapter  9,  for  example,  the  nearest  scattering  cosine  in  an  Ss  level 

symmetric  angular  quadrature  is  //£,  =  0.922  «  23°).  If  Oak  Ridge 

National  Laboratory’s  175  group  VITAMIN-J  group  structure  is  used  (see 
appendix  E),  109  of  the  175  groups  are  unable  to  downscatter  to  the  next 
lower  group  using  an  ordinate-to-ordinate  approach.  Lack  of  such  angular 
support  is  particularly  troublesome  for  within  group  and  next  group  scatter. 
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Mikols  provided  a  means  of  analyzing  such  behavior.  He  showed  where 
quadrature  and  group  structure  combinations  fail  for  neutral  particle 
scattering!!.  Using  the  same  type  of  analysis,  it  is  also  possible  to  show  that 
some  group-to-group  scatters  get  skipped  even  though  the  physics  says  they 
shouldn’t. 

The  difficulties  encountered  in  calculating  the  ordinate-to-ordinate  cross- 
sections  and  the  failure  of  the  method  for  sparse  angular  quadratures  have 
made  it  an  unacceptable  choice  in  the  transport  community. 


Non-negative  Functional  Scattering  Cross  Section  Representation 


Another  method  to  solve  the  source  negativity  problem  is  to  represent  the 
scattering  cross-section  with  a  strictly  non-negative  function.  One  such 
method  was  proposed  by  Dahl®. 

Dahl  also  attempted  to  calculate  the  group-to-group  scattering  cross 
section  as  a  function  of  angle,  and,  like  Odom,  Dahl’s  method  requires  storage 
of  the  angular  fluxes  during  the  transport  calculation.  Dahl  replaces  the 
conventional  finite  Legendre  scattering  cross  section  representation  with  an 
exponential  form 


f  L  ^  ^  ^ 

(Tg'g  (Q'-Q)  =  exp  ^  ^g'gPi  (O'  •  Q) 
W=o 


(37) 
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To  determine  each  X^g',g ,  the  constraint  that  the  original  cross  section 


moments,  a  ^g',g ,  be  preserved  is  imposed.  This  results  in  L+1  coupled  non¬ 
linear  equations  which  are  solved  for  the  X^g',g  ‘s: 


g'g 


(L  ^ 

U'=o 


=  0,  l  =  0,l,...,L. 


(38) 


Unfortunately,  Dahl  found  it  computationally  impossible  to  match  all  the 
moments  of  the  Legendre  polynomial  representation.  He  presented  an 
example  in  which  the  zeroth  through  fifth  moments  were  matched,  and 
showed  that  the  sixth  through  fifteenth  moments  of  the  exponential 
approximation  were  quite  close  to  the  corresponding  moments  of  the  actual 
cross  section.  While  this  is  promising,  he  has  not  presented  a  reliable  method 
for  finding  the  values  of  the  coefficients,  X’'g',g.  Dahl  warns  that  the  accuracy 
of  the  method  is  sensitive  to  the  number  of  equations  the  non-linear  solver 
must  evaluate.  If  too  many  input  moments  are  used,  the  exponential 
expansion  calculated  using  the  code,  XREP^^,  no  longer  preserves  the  original 
moments.  Thus,  there  is  an  inherent  limit  to  the  accuracy  of  the  method. 
Further,  while  the  method  is  strictly  positive,  it  should  only  be  non-negative. 
Even  where  there  are  scattering  transfers  which  are  strictly  not  allowed,  the 
exponential  representation  has  a  non-zero  value.  Finally,  even  if  this  method 
could  match  the  scattering  cross  section  exactly,  the  same  ordinate-to- 
ordinate  scattering  limitations  described  above  (lack  of  the  angular  support) 
apply. 
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Use  of  Piecewise  Constant  Functions 


The  methods  described  so  far  treat  the  scatter  as  being  only  from  one 
directional  ordinate  to  another.  Namely  corresponds  to  the 

scatter  from  ->  on  the  unit  sphere.  Borsari  suggests  an  alternate 
approaches,  jjg  suggests  thinking  of  the  transfer  as  being  the  average  scatter 
from  one  patch  of  solid  angle  on  the  sphere  to  another. 

Instead  of  the  directional  ordinates  representing  a  single  direction  in 
space,  consider  partitioning  the  possible  angular  directions  of  a  sphere  into 
(what  we  call)  facets.  With  each  facet  representing  a  separate  and  distinct 
patch  of  solid  angle  or  area  on  the  sphere.  Mathematically,  a  facet  on  the 
unit  sphere,  "U,  corresponds  to  a  portion,  AQ,  of  the  sphere’s  total  solid  angle. 
The  union  of  these  non-intersecting  patches,  or  facets,  cover  the  sphere,  so 
that 

iv 

^AQ„=jdQ  =  l.  (39) 

n=l  ‘ti 


The  scattering  from  ordinate-to-ordinate  (Q„.  Q^j)  is  now  represented  as 
an  average  of  the  scatter  from  facet  to  facet  ( AQJj-  AQ„  )^3_  And  the 
average  scattering  cross  section  is  defined  as 


(40) 
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To  accomplish  the  integration,  Borsari  assumed  a  piecewise  continuous 
function  (p.c.f.)  of  crggi/xj^),  so  that  <yg>g{ni)  is  approximated  as  stepwise 

bands  on  the  sphere  (Figure  4). 


Figure  4.  Borsari  scattering  cross-section  p.c.f. 

Borsari  spent  much  of  his  effort  looking  at  the  best  way  to  create  the  cross 
section  values  at  the  given  cr g gin l)  mesh  points  (referred  to  as  the  gamma 

grid).  He  examined  gamma  grids  that  depended  only  on  the  facet 
quadrature,  and  gamma  grids  that  depended  on  the  behavior  of  crggiML)- 

either  case,  the  accuracy  of  the  method  depended  on  the  resolution  of  the 
gamma  grid.  To  calculate  the  average  cross  section,  grid  values  were 
interpolated.  A  Cartesian  product  angular  quadrature  simplified  finding  the 
intersections  of  the  gamma  grid  with  facet  boundaries.  This,  in  turn,  made 
calculating  the  integral  of  equation  (40)  simpler,  but  restricts  Borsari’s 
method  to  Cartesian  product  quadratures. 
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The  implication  in  calculating  the  scattering  source  can  be  found  by 
returning  to  equation  (16) 


q^g(r,Q)  =  X  jdn'jdE\dE'cj^iE'  E,Q' -QWgiE')  (41) 


^-1 


where  the  source  is  now  approximated  as 


g'=ln'=l 

xj-^j-^jdEjdE'cT^(E'  E,Q'-Q.)Wg(E')  y/g'A^)- 


(42) 


The  scattering  source  term  takes  on  the  exact  same  form  as  equation  (19), 


G  N 


Qg,n(^)=  Z  'L^n’^g'gi^n'  '  ^n)¥ g' ,n' i^)  , 

g'=ln'=l 


(43) 


with 


and  the  angular  flux  unaltered. 

Not  discussed  by  Borsari,  but  mentioned  here  is  the  immediate  advantage 
of  this  technique  over  ordinate-to-ordinate  methods  in  that  it  guarantees 
angular  support. 

In  averaging  the  scatter  over  the  facets  some  of  the  scatters  will  take 
place  near  facet  boundaries.  As  and  approach  each  other  on  a  facet 
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boundary,  Thus,  from  equation  (34),  for  any  fixed  group  structure 

there  will  always  be  some  scatter  such  that 

Ojio.’  ^  q)  <  +  «l).  n'eQ^,,QeQ^.  (46) 

where  al  is  an  arbitrarily  small  number  greater  than  zero.  Such  behavior  is 
independent  of  the  facet  dimension  (or  quadrature  size)  Because  the  scatter 
occurs  at  the  facet  boundary.  Because  the  separation  between  neighboring 
facets  is  zero,  there  is  always  the  possibility  of  scattering  from  one  facet  to 
the  other  no  matter  how  small  the  loss  of  energy  in  the  collision.  And 
because  the  entire  sphere  is  tiled  with  facets,  every  group-to-group  scatter 
that  can  occur  (based  on  the  scattering  physics)  will  occur. 

The  Borsari  method  is  essentially  an  ordinate-to-ordinate  paradigm  that 
required  interpolation  schemes  over  his  values  when  scatter  spanned 
multiple  squares  on  his  Cartesian  grid.  The  intricate  relations  between  the 
way  the  source  facet  was  partitioned  for  integration  and  the  way  the  grid  of 
jUi  values  is  selected  affects  accuracy  and  non-negativity.  Solving  the 

integration  of  equation  (42),  Borsari  did  not  discuss  how  <T^gg{r,Hi)  is  found 

at  the  requisite  gamma  grid  points.  Further,  using  such  a  grid  required  that 
the  sphere  be  tiled  with  a  Cartesian  facet  quadrature.  This  made  calculation 
of  the  integrals  simpler,  but  the  method  can  not  modified  to  use  more 
symmetric  quadratures. 
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The  use  of  piecewise  continuous  functions  (or  facets)  also  required  the 
rebalance  of  the  cross  section  matrices  to  conserve  the  first  angular  moments 
of  the  scatter  (thus  preserving  the  currents).  Borsari’s  approach  requires 
that  the  rebalance  be  performed  during  execution  of  the  transport  code,  thus 
adding  another  level  of  complexity  to  the  method.  Further,  his  rebalance 
technique  may  destroy  so\u*ce  non-negativity ^3. 

While  Borsari’s  idea  of  the  facet  representation  holds  promise,  his 
implementation  is  complicated,  dependent  on  the  appropriateness  of  the  facet 
gridding  for  the  t5q)e  of  scatter,  and  wed  to  a  Cartesian  angular  quadrature 
type.  As  a  consequence,  most  of  the  paper  is  more  mathematical  rather  than 

practical.  No  effort  was  made  to  actually  calculate  for  real  materials. 

Present  methods  fail  to  adequately  address  the  problems  of  anisotropic 
scattering  in  neutral  particle  transport.  There  is  a  need  within  the  transport 
community  for  the  creation  of  a  new  approach  that  is  simple  to  calculate  and 
implement  in  transport  codes.  We  propose  a  completely  novel  approach  —  an 
approach  that  incorporates  the  sampling  of  the  physics  used  by  MCNP  at  its 
heart,  and  takes,  advantage  of  the  benefits  a  facet  method  has  over  ordinate- 
to-ordinate  techniques.  Our  goal,  as  always,  is  to  return  to  the  physics  and 
model  the  transfer  matrix  as  accurately  as  possible  without  resorting  to 
expansions  and  fix-ups. 
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The  Monte  Carlo  Facet  Method 


Return  to  the  physics  is  the  overarching  philosophy  of  this  approach. 

When  the  results  of  discrete  ordinates  codes  are  presented,  they  are  almost 
universally  benchmarked,  validated  and  otherwise  compared  to  similar 
problems  run  using  the  Monte  Carlo  transport  code  MCNP^'^.  MCNP 
stochastically  models  the  transport  process  as  continuous  in  space,  energy 
and  angle.  Experimentally  measured  cross  sections  are  combined  with 
predictions  from  nuclear  model  calculations.  The  results  of  this  combination 
are  incorporated  into  sets  of  data  such  as  the  Evaluated  Nuclear  Data  File 
(ENDF)  libraries^^.  MCNP  samples  this  data  directly  from  continuous  energy 
cross  section  information  obtained  from  the  ENDF  data.  The  continuous 
energy  cross  section  library  used  by  MCNP  is  known  as  an  ACE  library.  Of 
course,  discrete  ordinates  methods  are  discrete  in  space,  energy  and  angle. 
Our  goal  is  to  model  the  energy-angle  relationship  as  closely  to  the 
continuous  model  (much  like  MCNP)  as  is  possible.  To  this  end,  our 
technique  samples  directly  from  the  ACE  tables  in  much  the  same  way  as 
MCNP. 

Following  Borsari,  we  define  the  average  group -to-group  scattering  cross 
section  as 


dQ 

ACln 


jdE'jdE  a^iE'  E,Q'  •fi)W(E') 

g'  g _ 


\dE'W{E') 

g' 


(47) 
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where,  as  Borsari  did  ,  we  consider  the  average  scatter  of  particles  from  one 
facet,  AQ' ,  to  another,  AQ  as  shown  in  Figure  5. 


Figure  5.  Scattering  from  AQ'  AQ . 


Substituting  equation  (2), 


g" 


Unlike  Borsari,  we  chose  to  evaluate  the  integral  in  equation  (47)  by 
Monte  Carlo  integration.  Before  that  can  be  accomplished,  equation  (48) 
needs  to  be  put  into  a  form  that  is  readily  evaluated  using  MC  techniques. 

For  simplicity  the  numerator  and  denominator  of  equation  (48)  are 
examined  separately  and  re-written  as 


—Num 


1 


JdE'JdQ' 

g'  n' 


cr»(£;OW^(E')JdnJd£  (49) 

n  g 


and 
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(50) 


g' 

Recall  that  the  integration  over  all  scattered  energies  and  angles  of 

f{E'  ->  )  equals  the  expectation  value  of  the  number  of  neutrons 

that  emerge  from  the  collision  given  an  incident  energy,  E' .  The  same 

integration  over  the  scattered  group  and  facet,  namely 

^(Kl^dE  f(E'  E,Qn'  Clji)  is  the  expectation  value  of  the  number  of 
n  g 

neutrons  that  emerge  from  a  collision  and  scatter  into  energy  group  g  and 
facet  n. 

In  the  MC  facet  method  this  is  sampled  in  a  specific  way.  For  a  given 
material,  the  type  of  scatter  is  determined.  This  scattering  physics  might  be 
elastic,  inelastic,  thermal,  etc.  The  material  dependent  sampling  of  the 
physics  will  be  represented  with  the  variable,  r.  This  may  involve  sampling 
the  scattering  angle,  //(./«>  to  find  Hi,  and  E,  or  it  may  involve  sampling  E  to 
then  find  //£,  .  In  either  case,  based  on  the  physics  model  and  the  incident 
energy  the  lab  frame  scattering  angle,  and  scattered  energy,  E,  are 
determined.  The  lab  frame  azimuthal  scattering  angle,  co^,  is  sampled 
uniformly  in  2n .  The  scattered  direction,  Q ,  is  constructed  from 
Q',  and  coi .  Finally,  the  appropriate  scattered  group  and  facet  are 
determined  such  that  E  e  and  Q  e  AQ„ . 
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We  now  define  the  distributions  and  such  that 


cx)  1  7  2;r  7  00 


-00  -1 


U  0 


X  d{Mi.-/iL)d’‘(E-E)s[a-n)xi(B)x^(n) 


n  g 


(61) 


where  in  general 

00 

ldEg{E)s{E-E)^g(E), 

0 

JdQ^(Q)(52(Q_Qj  =  ^(Q), 

7i 

1 

J djui  gi/ui ) d{nL  -Ml)^  Sil^L ) . 

-1 

and  E,  Q,  and  constructed  by 

E  =  E{r,E'), 


and  Jii^=JijXt,E'). 

The  characteristic  functions,  and  Xn(^)  >  are  found  from 


and 


E  eAEg, 
else, 


Q  eAQ„, 
else . 


(52) 

(53) 

(54) 


(55) 

(56) 

(57) 


(58) 


(59) 


Substituting  equation  (51)  into  (49)  gives 
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(60) 


-zNum  1 


00  2rr 


JdE'JdQ’Jdr  (£)/?(«) 

g’  n'  -00 


Afi„.An„  ,  , 

^  ^  71  -00  0 


Applying  Monte  Carlo  sample  mean  integration  to  approximate  the 
integral  above  (Appendix  A),  the  integrals  in  t  and  coi  are  sampled  using 

the  distributions  shown.  The  integration  over  g'  is  found  by  sampling  E' 
uniformly  over  ^^Eg■  (represented  by  the  distribution  fg'{E')).  Similarly  the 

integral  over  n'  is  found  by  sampling  Q'  uniformly  over  AQ„. ,  or  . 

Showing  the  sampling  distributions  explicitly  in  (60)  we  have 


-Num  _ 


00  2;r 


^g'g,n'n 


AQ 


jdE'jdQ'Jdr  fdai,a'’mW(E')z§(E)z^(n) 

^  g'  n'  -00  0 

X  fA^')fn’{n')f,(T;E'.h')  , 


(61) 


and  for  the  denominator 


5^,.,,  =  AE^.jdE'  W(E')  4.(E') . 

g' 

The  resulting  Monte  Carlo  approximation  is 


(62) 


1  1 


—n'n  ^  i=l 

^g'g  ~  ,  M 


Mr  ^ 

'Z[aHEOW(.E;)z^(,Ei)Zn&i) 


'ZmEl)zi(.Ei)z^iai) 


M 


g,n  i=l 


M 


(63) 


where  M  is  the  total  number  of  incident  particles,  „  =  ^^Xg 

1=1 


the  number  of  scatters  to  group  g  and  facet  n,  and  the  sampling  of  incident 
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energy,  E'^ ,  is  done  coherently  for  the  two  summations  to  guarantee 
appropriate  normalization  over  the  selected  scatters.  The  sampling  is 
performed  as  follows: 

1.  E'i  is  sampled  uniformly  in  /^Eg , 

A 

2.  Q'-  is  sampled  uniformly  in  AQ„. , 

3.  El  and  Qj  are  found  consistently  from  the  scattering  physics  by 

sampling 

Typically  (for  elastic  and  level  inelastic  scatter)  the  sampling 
variable  r  represents  which  is  sampled  from  a  distribution  table 

based  on  the  selected  physics  model  and  E- .  From  ju^m 
Pi  is  calculated  using  scattering  kinematics,  and 
El  is  found  using  scattering  kinematics,  then 
0)1  is  sampled  from  fco^{o)i)  uniformly  in  [0,2;r), 
and 

Qj  is  constructed  geometrically  from  QJ ,  pi  =Qi  -Qi,  and  6)^ . 

4.  The  scattered  group  and  facet  are  rejected  unless  Ei  eAEg,  and 
Qj  e  . 

Implementation  of  the  MC  facet  method  is  simple  and  straight  forward. 
The  scattering  data  is  taken  directly  from  the  continuous  ACE  tables 
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constructed  directly  from  ENDF  data.  Any  energy  weighting  function  can  be 
used.  Algorithm  1,  below,  outlines  the  essential  steps  of  the  MC  facet 
method.  The  remaining  theory  sections  explain  how  each  step  is  performed. 
The  approach  is  completely  general.  To  demonstrate  it,  we  used  a 
particularly  convenient  angular  quadrature  and  implemented  only  a  few  of 
the  MCNP  scattering  laws. 


35 


Algorithm  1.  Basic  Algorithm  for  the  MC  Facet  Method 


Set  desired  facet  quadrature,  total  facets  =  N. 

Set  desired  group  structure,  total  groups  =  G. 

Do  g'=l..G 

Do  n'=l..N 

TallyNum  =  0  Array  Dimension  (NxG) 

TallyDen  =  0  Array  Dimension  (N  x  G) 

TallyErr  =  0  Array  Dimension  {N  x  G) 

Do  i  =  1..M  (until  converged  or  some  max  sample  size  reached) 
Draw  E'  uniformly  in  /sEg . 

DrawQ'  uniformly  in  . 

Draw  E  and  consistently  from  scattering  physics. 
-(Incorporates  MCNP/ACE  models  of  ENDF-B  data) 
Draw  (Oi^  uniformly  in  [0,2;::). 

Construct  Q  geometrically  from  Q' ,  //£,,  and 
Find  g  such  that  E  e  . 

Find  n  such  that  Q  €  AQ„ . 

TallyNum(n,g)  =  TallyNum  (n,g)  +  cj^{E')W{E') 
TallyDen(n.,g)  =  TallyDen(n,g)  +  W[E') 

TallyErr(n,g)  =  TallyErr(n,g)  +  \cr\E')W{E''^ 

End  Do 

,  M 

^g'g  «  AQ„.  [TallyNum(n,^)/TallyDen(n,g)] 


Err”.^  «1.96AQ„. 


Save  <jgg  and  Errgg  to  file  for  n  =  l..iV,  andg  =  1..G. 

End  Do  (Incident  Facets) 

End  Do  (Incident'Groups) 


1 

r  1  1 

r  1  n 

2] 

— TallyErr(7i,g) 

- 

— TaUyNum(n.,g) 

-^|^^TallyDen(n,g) 

§1^ 
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Advantages  of  the  MC  Facet  Method 


The  Monte  Carlo  facet  method  has  a  number  of  advantages  over 
previously  discussed  methods: 

Guaranteed  Non-Negative.  Equation  (48)  is  approximated  using  a  Monte 
Carlo  method,  where  the  values  are  accumulations  of  non-negative  values. 
The  scattering  matrix  is  zero  where  scattering  is  not  allowed  and  positive  for 
values  where  scattering  is  allowed. 

Guaranteed  Within  Group  Scatter.  Adjacent  facets  will  always  scatter  into 
each  other  no  matter  how  small  the  energy  loss  (equation  (46)).  This  means 
that  particles  can’t  be  trapped  in  a  group  (as  in  the  ordinate-to-ordinate 
method)  because  of  the  group  structure  used  no  matter  how  sparse  the 
angular  quadrature.  Similarly,  down  scatter  groups  can’t  be  skipped  as  an 
artifact  of  the  angular  quadrature.  This  will  be  demonstrated  in  chapter  9. 

Adjustable  accuracy.  Unlike  the  method  proposed  by  Dahl,  which  fails  if 
one  attempts  to  match  too  many  moments,  the  user  can  achieve  the  desired 
cross-section  accuracy  by  simply  increasing  the  number  of  trials  in  the 
calculation. 

Arbitrary  Quadratures.  Because  particles  are  tracked  directly  from  one 
facet  to  another,  any  shaped  facet  or  quadrature  structure  may  be  used.  This 
is  distinct  from  Borsari’s  analysis  which  applied  only  to  Cartesian  type  grids. 
Hence,  the  quadrature  can  be  triangular,  icosahedral,  etc. 
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Incorporates  ENDF  Scattering  Laws.  A  number  of  techniques  (e.g., 
Odom)  consider  only  elastic  and  level  inelastic  scattering  because  the 
integration  of  the  scattering  kernel  is  simplified.  Other  methods,  such  as 
Borsari,  make  no  attempt  to  address  the  type  of  scattering  involved.  In  fact, 
such  scattering  can  he  quite  complicated.  There  are  a  number  of  different 
ways  to  describe  the  t5q)es  of  scattering  involved.  The  Monte  Carlo  facet 
method  is  able  to  take  advantage  of  all  the  ENDF  scattering  laws  in  exactly 
the  same  way  that  MCNP  does.  Any  neutral  particle  scattering  law 
(neutrons  or  photons)  can  follow  the  exact  same  formalism  developed  here. 
All  the  scattering  possibilities  are  accounted  for  in  one  integration  of  the 
average  cross  section  (unlike  Odom’s,  for  example,  which  requires  a  separate 
integral  calculation  for  each  type  of  scatter).  Taken  together,  these 
advantages  suggest  a  significant  breakthrough  in  the  attempt  to  account  for 
anisotropic  scatter  in  neutral  particle  transport. 
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5.  Implications  of  Facet  Transfers  on  Neutral  Particle  Transport 
This  chapter  outlines  some  mathematical  cautions  concerning  application 
of  the  average  scattering  technique  to  the  transport  code.  It  ends  by 
demonstrating  how  this  method  is  used  with  TETRAN,  which  is  the  neutral 
particle  transport  code  developed  at  AFIT^. 

Effect  of  Scatter  Averaging 

Recall  that  in  discrete  ordinates  theory  the  scattered  angular  flux  is 
approximated  by 

G  N 

■Q,^)Wg,(E')y/gr^,(f)  (64) 

g'=ln-=l  ^  ^ 

where  ^g'^n>  is  the  total  group  angular  flux  evaluated  at  a  specific  directional 

ordinate.  Here,  evaluation  of  the  group-to-group  scattering  source  represents 
the  explicit  transfer  of  particles  from  all  incident  groups  to  one  scattered 
group  and  all  incident  directions  to  a  specific  directional  ordinate. 

Use  of  the  average  scattering  cross  section,  cF”-” ,  puts  a  small  but 

important  conceptual  twist  on  this  transfer.  Here  the  scattering  from  all 
incident  directions  to  a  specific  scattered  directional  ordinate  is  approximated 
as  the  average  scatter  of  all  incident  directions  to  a  facet  surrounding  (in  an 
appropriate  way)  the  scattered  directional  ordinate. 


39 


One  advantage  of  such  a  method,  previously  discussed,  was  that  this 
guarantees  that  all  group-to-group  scatters  will  have  somewhere  to  scatter — 
regardless  of  the  scattering  anisotropy.  The  disadvantage  of  such  a  method  is 
that  it  can  contribute  to  angular  dispersion  of  the  scatter  and  that  it  fails  to 
conserve  the  first  angular  moment  of  the  scatter,  which  is  important  for 
diffusion  problems.  We  examine  the  dispersion  effects  in  subsequent 
chapters  and  discover  that  their  impact  on  the  transport  results  is  minimal. 

In  the  following  chapter  we  address  how  to  redistribute  the  scattering  matrix 
in  order  to  conserve  the  first  angular  moment  of  the  scatter. 


Standard  MC  Facet  Method  Implementation 
To  speed  the  calculation  of  the  scattering  source 


G  N  ^ 

~  S  g' a' 

g'=ln'=l 

we  define  the  group-to-group,  facet-to-facet  transfer  matrix,  T,  as 


(65) 


rpnn 


=  W 


^Vg  =  aq,. 


-z/i'n 
^g'g  ■ 


(66) 


Although  T  is  an  array  of  rank  4,  we  think  of  it  as  a  matrix  because  only 
the  rank  two  array  section  ,  for  one  combination  of  groups,  ^  and  g,  need 
be  loaded  into  memory  during  the  transport  code  run  at  any  given  time.  This 
allows  the  transport  code  to  perform  a  direct  matrix  multiply  using  an 
optimized  intrinsic  routine  in  Fortran-90.  The  scattering  source  becomes 
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(67) 


G  N 
g'==l  n'=l 


The  transport  code,  TETRAN,  used  for  this  research  only  allows  for 
downscatter,  so  equation  (67)  is 


(68) 

g'=ln'=l 


Using  the  equations  for  and  cr”-”  from  the  previous  chapter,  and  because 
AQ„.  =  AQ„  for  the  Cartesian  quadrature  we’ve  selected,  the  T-matrix  is 


rpn'n  ^  gt 


^  -  1 

2k“(£;)»'(E,o4(Ei)/?(5i) 
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M 
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Finally,  for  ease  of  notation  when  describing  error  calculations,  we  define 


Mr  2  -  ■ 

2  xf(.E0x2&i) 

i=iL 

^g'g  ““ 


M 


M 


(70) 


EH'(E04(E;);ir?(Si) 

1  =  1 


In  later  chapters,  we  examine  the  behavior  of  the  T-matrix,  and 
demonstrate  its  application  in  a  number  of  transport  calculations. 
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6.  Description  of  the  Monte  Carlo  Transport  Approach 


Use  of  Monte  Carlo  (MC)  integration  is  advantageous  for  two  reasons. 
First,  Monte  Carlo  methods  tend  to  be  more  efficient  for  integrals  of  greater 

A 

than  a  few  dimensions  (here  there  are  five:  E',  and  cdi  ). 

Second,  the  way  the  MC  draws  are  performed  to  calculate  the  scattering  cross 
sections,  angular  distributions  and  neutron  energies  after  collision  follows 
the  physics  of  the  scatter  and  closely  resembles  the  methodology  of  Los 
Alamos  National  Laboratory’s  (LANL)  Monte  Carlo  transport  code  MCNP. 

The  previous  chapter  briefly  described  how  is  calculated.  This 

chapter  details  the  implementation  of  the  steps  in  our  code,  T-Scat.  T-Scat 
models  the  physics  of  the  scattering  kernel  as  described  in  the  previous 
chapter  and  performs  the  requisite  Monte  Carlo  calculations.  It  is  based  on 
MCNPi4  and,  as  such,  its  capabilities  can  be  greatly  expanded.  This  research 
examined  only  the  scatter  of  neutrons,  although  the  same  formalism  holds  for 
gammas  as  well. 

Where  variables  are  described  below,  it  can  be  assumed  that  they  come 
directly  from  the  ACE  cross  section  library.  Where  exact  tabular  values  are 
not  available,  linear  interpolation*  is  typically  applied. 


*  Some  tables  (typically  of  larger  nuclides)  require  other  types  of  interpolation  (e.g.,  Log- 
Log).  Such  interpolation  is  not  required  for  the  hght  nuclei  examined  here. 
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Sample  the  Incident  Energy 


Once  the  ACE  library  is  read  in,  the  first  step  is  to  draw  E'  uniformly 
from  the  incident  energy  group.  First  determine  E'froin 

E'  =  E'g^^+p{E'g-E’g^-^)  (71) 

where  /?is  a  random  number  on  the  interval  [0,1). 

Sample  Incident  Direction 

Facet  Quadrature 

A  preliminary  step  is  to  define  the  tiling  of  the  sphere,  “U,  into  facets,  AQ„  . 

This  defines  the  angular  quadrature.  This  effort  used  a  modified  Cartesian 
grid  with  polecaps.  While  any  facet  quadrature  (or  surface  tiling)  can  be  used 
(an  icosahedron  for  example)  the  Cartesian  grid  makes  sampling  of  the 
incident  facet  simple  and  determination  of  the  scattered  facet  straight 
forward.  A  74  facet  quadrature  is  shown  in  Figure  6  and  Figure  7. 


.  Figure  6.  74  Facet  Polecap  Quadrature 
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Figure  7.  74  Facet  Structure  on  a  Cartesian  map. 


Most  transport  codes  are  unable  to  use  this  quadrature  set  because  there 
are  ordinates  at  the  poles.  The  transport  code  used  for  this  research, 
TETRAN,  uses  an  unstructured  tetrahedral  mesh  and  can  use  any  ordinates. 
The  polecap  quadrature  facilitates  rebalancing  to  conserve  the  first  angular 
moment  of  the  scatter,  as  described  in  the  next  chapter.  To  be  general, 
T-Scat  can  also  create  scattering  matrices  for  Cartesian  meshes  without 
ordinates  at  the  poles,  but  no  rebalancing  technique  has  been  developed  for 
such  a  quadrature  set. 

Each  facet  of  the  sphere  has  the  same  solid  angle.  The  azimuth  is  divided 
with  equally  spaced  lines  of  longitude  (A<y=  constant).  If  is  the  number 

of  longitudinal  lines  dividing  the  azimuth  then 
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(72) 


Act)  = 


‘In 


N. 


(O 


The  same  equi-partition  is  applied  to  the  polar  angle  with  a  different 
treatment  at  the  poles.  To  assist  in  the  conservation  of  the  first  angular 
moment  of  the  scatter,  a  cap  has  been  placed  on  each  pole.  This  polecap  has 
the  same  solid  angle  as  the  other  rectangular  (in  co  coordinates)  facets.  The 
result  is  that  the  polecap  has  a  smaller  polar  width  than  the  other  facets. 

Dividing  the  polar  angle  [-1,1]  with  divisions  the  polecap  polar  width 


IS 


2  +  JV„(Nj-2)’ 


(73) 


while  for  a  rectangular  facet 


2  +  iV^(iV^-2) 


(74) 


The  total  number  of  facets,  N,  is 


N  =  2  +  N^(N^-2). 


(75) 


As  a  means  of  simplifying  the  description  of  the  quadrature,  this  paper  refers 
to  a  given  quadrature  set  as  x  ->  N  so  that  8  x  12-»74  defines  a 


45 


quadrature  set  that  has  6  x  12  =  72  rectangular  facets  and  two  polecaps  for  a 
total  of  74  facets  or  angular  ordinates. 

In  general,  the  directional  ordinates  for  any  arbitrary  facet  quadrature 
are  found  from 


f  h 

f  dQQ 

_  Jn 

dCl  ^ 

Ilf  dClQ 

Wm 

where  the  directional  vector,  Q ,  is 

cos{co)i  +  sin{o))  j  +  , 


(76) 


(77) 


and 


!q 


is  the  Euclidean  norm  of  Q , 


=  ./Q?+Qf +q| 


(78) 


Integrating  equation  (76)  over  a  given  facet  where  co  e  \a)_  ,co+  ] ,  and 
4e[^-,^+]  gives 


p  r -  ^  - 

I  (iQQ=  cos(©)i sin(a))j  +  ^k 

6). 


(79) 


The  azimuthal  angle,  relative  to  the  facet  edge,  is  easily  found  as 


0)  =  tan 


-1 


Jdfi?sin(«)/  Jdfi?cos(ffl) 


=  (0_+-{(0^-(0_). 


(80) 
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The  polar  cosine  is  a  bit  more  difficult  to  calculate.  The  k  component  of 
equation  (76)  is 


(81) 

+  4  sin  2  ](^+  +  sin'^  )  -  sin"^  (^_  )j 


For  even  modest  facet  widths,  Aco  <nl  A,  this  is  approximated  by 


+  ^(^+)-sin 


1/2  (82) 


and  for  <  0.3,  this  is  further  simplified  by 


iJ±±lA-P  M 

2  2  ’ 


(83) 


which  is  good  to  about  one  percent.  For  the  8  x  12-^74  facet  quadrature  set, 
the  polar  angles  of  the  ordinates  using  equation  (81)  are  found  at 
4^  =  ±  0.165, ±  0.494,±  0.826,±  1.00 .  While  the  polar  angles  using  equation  (83) 
are  found  at  ^  =  ±  0.162, ±  0.486,±  0.811,±  1.00 .  For  this  research, 

equations  (80)  and  (81)  are  used  to  obtain  directional  ordinates. 

There  is  a  substantial  literature  discussing  more  complicated  two-angle 
quadrature  sets  particularly  in  modeling  molecular  structures  Some 

of  these  could  be  drawn  from  to  create  more  appropriate  quadrature  sets. 
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The  Monte  Carlo  facet  method  could  be  used  for  such  quadratures  as  shown 
in  Figure  8. 


(C) 


Figure  8.  Quadrature  schemes  proposed  by  Lebedev^®,  a)  Lebedev 
quadrature,  b)  Uniform  triangulation,  c)  Chebyshev  quadrature 

Facet  Sampling 

The  incident  facet  is  sampled  in  the  following  way.  Let  pj  and  P2  again  be 
random  numbers  on  [0,1).  Then  for  the  given  incident  facet,  n', 


(84) 

and 

0)'  =  0>'_+  p2(.G>+  -0)-). 

(86) 

These  values  are  stored  until  the  lab  frame  scattering  angle  is 
determined. 
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Select  the  Collision  Nuclide 


There  are  two  ways  multiple  nuclei  (or  isotopes)  in  a  material  may  be 
handled.  In  the  first  method,  the  scattering  matrices  for  the  nuclides  (or 
isotopes)  are  calculated  separately.  The  resultant  matrices  are  then  mixed  in 
an  appropriate  way  after  the  calculation.  In  general,  when  there  are  M 
different  nuclides  forming  the  material  then 

M 

Tu^=Y.Tih,  (86) 

1=1 

where  /j  is  the  atom  fraction  of  nuclide  i.  For  H2O  this  is 

'^H^O  =^Th+Tq.  (87) 

The  alternative  is  to  have  T-Scat  select  the  collision  nuclide  as  part  of  the 
sampling  of  the  physics.  When  this  is  done,  T-Scat  follows  the  methodology 
of  MCNP.  If  there  are  M  different  nuclides  forming  the  material  then  the 
nuclide  is  selected  if 

k-1  M  k 

fi  <  fi  ^  Z^i  fi  > 

i=l  1=1  i=l 

where  ^  is  a  random  number  on  [0,1),  and  ct\  is  the  total  cross  section  of 
nuclide  i.  The  total  and  scattering  cross  sections  are  found  as 

M 

<T‘(E')  =  £o-‘(E')/i.  (89) 
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and, 


(90) 


M 
i=l 

While  the  second  method  provides  the  user  with  a  scattering  matrix  that 
doesn’t  need  further  mixing,  the  first  method  has  a  number  of  advantages.  If 
five  million  draws  were  used  to  create  each  of  and  then  by 

mixing  Tjj  with  To  the  number  of  samples  has  effectively  doubled.  Further, 

once  the  scattering  matrices  have  been  created  for  a  number  of  selected 
isotopes,  they  form  a  fixed  data  library.  The  user  can  then  select  desired 
nuclides  to  build  the  material  of  choice.  This  is  comparable  to  present  cross 
section  libraries  ENDF/B-VI  and  MATXS-10  discussed  in  detail  later. 

Determine  Type  of  Scatter 

For  a  known  nuclide  and  incident  neutron  energy  in  the  lab  frame,  E',  the 
next  step  is  to  determine  if  the  scatter  is  elastic  or  inelastic.  We  proceed  in 
the  following  way:  Selection  of  an  elastic  collision  is  made  with  probability 

<7,,.(£;')  +  <Trf(E’)  <TT(E')-a^(E-) 

where 

a^i  is  the  elastic  scattering  cross  section, 

<Jl^  is  the  inelastic  scattering  cross  section, 

cTq  is  the  absorption  cross  section,  and 

(Tt’  is  the  total  cross  section,  cry  =  cTgi  +  cxij^  +  cr^ . 
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For  an  incident  neutron  energy  (in  the  lab  frame),  E',  bounded  by  two 


values  on  the  cross  section  table  E'k<E'<E'k+l,  where  k  is  the  table  index,  the 
cross  section,  cy(E'X  is  found  from 

<7(£')  =  )-a(E'k  )] .  (92) 

E  k+\~E  k 

A  plot  of  the  total  collision  scattering  cross  section,  a'f{Ej^),  taken  from  the 

continuous  energy  ACE  file  for  oxygen  is  shown  in  Figure  9  and  Figure  10. 
The  plot  of  Figure  10  shows  the  detail  of  the  resonance  region  for  oxygen. 


Total  Collision  Cross  Section 


Energy  (MeV) 


Figure  9.  Total  collision  cross  section  from  ACE  file. 


Total  '®0  Collision  Cross  Section 


Figure  10.  Total  collision  cross  section  from  ACE  file  (1-20  MeV) 
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If  the  collision  is  determined  to  be  inelastic,  the  type  of  inelastic  reaction, 
n,  is  sampled  from 

(93) 

1=1  i=l  i=l 

where  N  is  the  number  of  inelastic  reactions  [e.g.,  (n,2n),  (n,nl),  (n,n2),  etc.] 
with  cross  sections  found  from  the  ai ’s.  The  same  random  number  used  to 

determine  if  the  scatter  was  elastic  is  now  re-normalized  to  the  interval  [0,1) 
to  select  the  type  of  inelastic  collision.  This  new  random  number,  p,  is  found 
from  the  old  random  number  with 


Pold  -  Pel 
l-Pei 


(94) 


thus  saving  another  call  to  the  random  number  generator. 

Just  as  in  MCNP,  for  elastic  and  discrete  inelastic  scatters,  the  scattered 
angle  is  determined  from  equi-probable  angular  distribution  tables,  and  the 
scattered  energy  is  determined  from  two-body  kinematics.  For  other  inelastic 
processes,  the  energy  of  the  scattered  particle  is  determined  from  secondary 
energy  distribution  laws,  which  vary  according  to  the  particular  inelastic 
collision  modeled^^.  To  simplify  development  of  the  code  while  still 
demonstrating  the  potential  of  the  MC  facet  method,  T-Scat  limits  its 
calculations  to  elastic,  discrete  inelastic,  and  continuum  scattering  using  the 
inelastic  evaporative  spectrum  scattering  law  (described  below).  The  other 
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ENDF  scattering  laws  used  by  MCNP  can  be  easily  added  to  the  T-Scat  code 
as  simple  case  options  in  FORTRAN-90. 

For  both  elastic  and  level  inelastic  scattering  t  is  the  scattering  angle  in 
the  center-of-mass  frame,  /icm  •  Thus,  ju^n  is  sampled  from  its  tabulated 
distribution,  and  it  determines  jui  and  E.  The  ju^n  angular  distribution 
tables  are  represented  in  the  same  way.  The  tables  consist  of  32  equi- 
probable  cosine  bins  in  the  center-of-mass  frame  given  at  a 

number  of  different  incident  lab  frame  energies,  E’j ,  where  j  is  the  table 
index.  (Typically  the  table  for  the  lowest  incident  energy,  E'q  ,  for  the 

selected  reaction  is  assumed  isotropic.)  The  table  format  is  outlined  in  Table 
1.  Note  that  there  are  33  scattering  angles  defining  the  32  bins. 


Table  1.  Equi-probable  scattering  bins  from  ACE  file. 


Table 

Energies 

CM  Scattering  angles 

E’q 

... 

//""o,32 

* 

* , 

E'j 

... 

32 

The  table  values  for  hydrogen  are  shown  in  Figure  11,  while  the  same 
tjq)e  of  chart  is  shown  for  '^Li  in  Figure  12. 
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Hem  Equi-Probable  Bins 
Elastic  Scatter 


Figure  11.  Equi-Probable  Bins  Elastic  Scattering  Table 


Energy  (MeV) 


Figure  12.  "^Li  Equi-Probable  Bins  Elastic  Scattering  Table 


Note  that  the  table  for  hydrogen  is  nearly  isotropic  for  all  energies  in  the 
center-of-mass  frame.  Hydrogen  gets  its  anisotropy  when  the  scatter  is 
converted  to  the  lab  frame  because  of  its  low  atomic  weight.  The  anisotropy 
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in  the  center-of-mass  frame  is  seen  at  high  energy  for  lithium.  The  data  table 
can  be  used  to  create  a  cumulative  distribution  function  of  the  scatter  at  the 
given  energies.  Such  a  table  is  shown  in  Figure  13  for  ^Li.  From  the  figure, 
it  is  easy  to  see  that  at  high  energies  the  probability  of  drawing  a  scatter  in 
the  forward  direction  (//cm  =  1)  is  greater  than  at  low  energies.  This  forward 
distribution  is  enhanced  when  converting  the  scatter  to  the  lab  frame  for 
light  nuclides. 


If  E'j  is  the  energy  of  table  j  and  E'j^i  is  the  energy  of  table  j+1  then  a 
value  of  the  scattering  angle  //„„  is  sampled  from  table  j+1  with  probability 
-Ej'j  and  fi*om  table  j  with  probability  -E'jY 

A  random  number  p  on  the  interval  [0,1)  is  then  used  to  select  the  cosine 
bin  such  that  i=floor[32p  +1].  The  value  of  //cm  is  then  computed  as 

Mem  =  Mi  +  (32p - iXPi+i  -Mi)-  (95) 
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If  is  isotropic  it  is  found  on  [-1,1)  as 

f^cm  ~  ~  1  • 


(96) 


Elastic  Scattering 

For  elastic  scattering,  once  the  scattering  angle  is  calculated  the 
scattered  energy  E  is  determined  from  two-body  kinematics.  In  the 
laboratory  frame  it  is 


E  =  E' 


+2Ai^  cm 

(1  +  A)2 


(97) 


where  A  is  the  atomic  weight  and  E’  is  the  incident  neutron  energy. 


Level  Inelastic  scattering 

For  level  inelastic  scattering  (a.k.a.  MCNP  law  3),  the  energy  of  the 
scattered  neutron  is  found  from  basic  kinematics.  In  the  center-of  mass 
frame  the  scattered  energy  is 


^cm  - 


A  Y 

L,  l<?M+i)l 

U+iJ 

A 

(98) 


where  Q  is  the  excitation  energy  of  the  nucleus.  This  equation  can  also  be 
used  to  find  the  exiting  energy  in  the  center  of  mass  frame  for  elastic 
scattering  by  setting  Q=0.  The  outgoing  energy  in  the  laboratory  frame  is 

E'  +  2jUc,n(A  +  l)4^, 


E  =  Ecm+- 


'cm 


(A+iy 


(99) 
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Evaporation  Spectrum 

The  evaporation  spectrum  reaction  is  one  way  to  model  scattering  to  the 
continuum.  The  distribution  is  sampled  independently  of  the  angular  draw. 
Hence,  energy  conservation  is  not  strictly  enforced  for  each  individual 
collision,  but  overall  energy  conservation  is  approached  as  the  number  of 
draws  becomes  large.  The  spectrum  is  given  as^^ 

f{E'  -^E)  =  ,  (100) 

where  the  nuclear  temperature,  T(E'),  is  a  tabulated  function  of  the  incident 
energy.  The  energy  U  is  provided  in  the  library  and  is  assigned  so  that 
is  limited  by  0  <  <  E'  -U  .  The  normalization  constant  C  is  given  by 

C"^  =  T2[i  -  (1  +  {E'  - 1/)  /  T)] .  (101) 

As  in  MCNP,  T-Scat  samples  the  density  function  with 

E„„=-r(E')ln(/?iP2),  (102) 

where  pi  and  p2  are  random  numbers  on  (0,1]  and  are  rejected  if 
Ecm  ^  E’  —  XJ . 

This  completes  the  set  of  scatter  reactions  that  are  currently  implemented 
in  T-Scat.  There  are  some  thirteen  different  laws  for  neutrons  and  more  for 
photons  (See  Appendix  B;  The  Inelastic  Neutron  Scattering  Laws  of  MCNP). 
All  could  be  incorporated  into  the  same  algorithm  with  the  goal  of 
determining  the  scattered  angle,  pi ,  and  scattered  energy  E. 
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Finding  the  New  Direction 


Once  the  scattering  angle  (in  the  center-of-mass  frame)  and  energy  (in 
both  center-of-mass  frame  and  lab  frame)  are  known,  the  direction  of  the 
exiting  neutron  is  found  in  the  following  way  where  we  assume  that  the 
target  is  at  rest.  First  the  cosine  of  the  scattering  angle  in  the  laboratory 
frame,  ,  is  found  from 


E. 


cm^ 


+  ■ 


E  A-flV£: 


(103) 


The  incident  particle  direction  cosines,  (ju',  t]',  ^')  are  rotated  to  new 
scattered  cosines  ( fx,  rj,  ^),  through  a  polar  angle  whose  cosine  is  and 
through  an  azimuthal  angle  sampled  uniformly.  Conceptually,  we  sample 
(Oj^  as  discussed  above.  However,  to  avoid  expensive  evaluations  of 
trigonometric  functions,  we  sample  its  components  by  rejection  as  follows. 

For  random  numbers  p-^  and  on  the  interval  [-1,1)  with  rejection  criterion 


(y  cy 

/?!  -I-  P2  - 1  >  the  general  equations  are 


p  =  p'pi^ 


V  =  Vi'liL  + 


Pl{PiV' -  P2P') 


(104) 


(105) 


and 


^  =  ^'pl  + 


■Pl) 

(pi  +pV\ 

1 

(106) 
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If  l-^'2  »  0,  then  an  alternative  approach  is  needed  to  avoid  catastrophic 
cancellation: 


and 


I  =  n'i 


#=^>L+- 


uliPlP'V'  +  P2^') 

Plj 

L 

/(pl +P2) 

^l(^-  -  P2M') 

(Pl 


(107) 


(108) 


(109) 


If  the  scattering  distribution  is  isotropic  in  the  lab  frame,  it  is  possible  to 
use  an  even  simpler  formulation  that  takes  advantage  of  the  exiting  direction 
cosines  being  independent  of  the  incident  direction  cosines.  In  this  case, 


P  =  2(/7i  +pI)-1, 

(110) 

1 1-p^ 

22’ 

(111) 

VPl  +P2 

.  |1-P^ 

^  =  2  2- 

(112) 

VPl  +P2 

where  and  are  rejected  ii  p\^r  p\>l. 

A  simple  variance  reduction  technique  would  be  to  split  the  scattered 
neutron  for  a  fixed  number  of  azimuthal  angles,  rather  than  draw  the 
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azimuthal  angle  for  a  single  number.  For  example,  select  a  random  number, 
p,  on  [0,1).  Then,  for  i=l,2,..360  select  360  azimuthal  angles  by 

+  .  (113) 

‘  ^  ^  360 

Now  that  the  physics  has  been  applied  to  the  scattering,  the  next  step  is  to 
identify  the  facet  and  group  in  which  Q ,  and  E  lie. 


Identifying  Scattered  Energy  Group  and  Facet 

As  required  by  equation  (63),  for  each  sample,  once  the  scattered  energy 
and  direction  are  known  the  appropriate  scattered  group  and  facet  must  be 
identified.  In  many  cases,  particularly  at  high  energies  for  heavy  nuclei, 
little  energy  is  lost  in  the  collision.  An  example  of  this  for  oxygen  is  shown  in 
Figure  14.  This  suggests  that  using  the  previous  scattered  energy  as  a  first 
guess  for  the  next  search  might  improve  performance. 


Figure  14.  Scattered  energy  spectrum  for  for  incident  E'  e  [17.333, 19.64]. 
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Once  the  scattered  energy,  E,  is  known,  a  search  is  performed  on  the 


group  structure  to  determine  the  appropriate  scattered  group 

El  6  (Eg,Eg_i].  A  bisection  search  can  be  done,  but  we’ve  found  that  at  high 

energies  using  the  last  known  scattered  energy  as  a  first  guess  improves  code 
performance  by  about  2%  (175  group  structure).  The  hunt  routine  was  found 
in  Numerical  Recipes^o. 

The  next  step  is  to  find  the  facet  for  which  Qj  e  AQ,j  .  For  the  Cartesian 
grid,  no  search  really  needs  to  be  performed.  Once  the  scattered  direction 
Qj  =  {jui,j}i,^i)  is  determined,  the  azimuthal  angle,  co,  is  found  from 


(Oi  =  tan 


-1 


N 


,  0)1  e 


[-7t,7t]. 


(114) 


The  azimuthal  facet  index,  ,  is 


ico  =  inti 


N. 


(O 


0)+  K 

2n  . 


+  1, 


while  the  latitudinal  facet  index,  is 


=  inti 


-f-  +  0.5iV,  1  +  1, 


The  scattered  facet  index,  n ,  is  then 


n  =  i^  + 


(115) 


(116) 


(117) 


For  the  given  incident  group,  g'  ,  and  facet,  n' ,  the  scattered  group,  g, 

and  facet,  n,  have  now  been  determined.  The  following  tallies  are 
accumulated: 
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(118) 


’—Nuin 


(new)  =  (old) + W'(£' ) , 


’—Den 


and  (^g^^’ni^ew)  =  erf J^,„(o/rf)  +  [cT(£:f)Ty(£f)]' 


After  a  set  number  of  draws,  M,  the  T-matrix  is  calculated  from 


(119) 

(120) 


rpn'n 

^g’g 


ACln'^gfl  (^gj'i'njnew) 
AQ„  M  ’ 


(121) 


where  Mg^  is  the  number  of  counts  in  group  g  and  facet  n.  The  T-matrix 


error  is  found  from 


n'n 

^g'g 


«L96 


V 


r~( 

' 

2] 

1  —  Err  (  \ 

-^(^g'gMnew) 

- 

1  -=Num 
^g'g/i'n 

> 

AQ. 


M 


gA 


(122) 


If  the  desired  error  is  not  met,  or  the  maximum  number  of  draws  has  not 
been  exceeded,  the  number  of  draws  is  doubled  and  the  process  is  repeated. 

A  number  of  data  files  are  output  including  the  T-matrix,  the  error  in  the 
T-matrix,  the  counts,  Mg^ ,  for  each  incident  to  scattered  facet  and  group,  as 

well  as  the  average  angle  of  scatter,  Qg'g ,  which  is  used  to  conserve  the  first 
angular  moment  of  the  scatter. 
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7.  Conserving  the  First  Angular  Moment  of  the  T-Matrix 


In  the  Monte  Carlo  facet  method,  scattered  neutrons  are  tallied  according 
to  the  facet  and  energy  group  into  which  they  scatter.  In  the  transport  code, 
the  angular  flux,  y/ g  „ ,  is  discretized  such  that  the  neutrons  propagate  along 

the  directional  ordinate,  .  When  calculating  the  average  scatter  in  the 

MC  facet  method,  in  the  case  of  highly  anisotropic  scattering,  only  a  small 
part  of  a  scattered  facet  may  be  reached.  The  assumption  that  these 
neutrons  travel  in  the  direction  of  that  facet’s  ordinate  (equation  (76))  is 
tantamount  to  a  very  large  redistribution  of  the  scattered  neutrons  away 
from  their  intended  direction. 

Here  we  fix  this  redistribution,  to  first  order,  by  actually  tracking  the  first 
angular  moment  of  the  scatter  (or  average  scattering  angle)  within  each 
angular  facet.  As  is  shown  below,  conserving  first  angular  moments  is 
required  if  the  T-matrix  is  to  support  codes  running  diffusion-like  problems. 

In  this  chapter  we  describe  the  diffusion  theory’s  requirement  for  first 
angular  moment  conservation.  We  then  describe  the  theory  behind  the 
rebalance  of  the  T-matrix.  Impact  of  this  rebalance  technique  on  transport 
computations  is  shown  in  subsequent  chapters. 

The  Pi  Approximation  and  Diffusion  Theory 

This  discussion  follows  the  presentations  of  Bell  and  Glasstone^,  and 
Marchuck  and  Lebedev^i.  It  is  a  brief  review  of  diffusion  theory,  and  helps 
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explain  why  conserving  the  average  angle  of  scatter  is  important  in  diffusion 
theory. 

The  steady  state  Pi  equations  used  in  diffusion  theory  are 
V  •  J(r ,  E)  +  C7(r ,  E)(j){r,  E)  =  (r,  E)  +  jdE'  ctsq  if,E'^E)<p(r,E'),  (123) 

and 

i V^(r , jE)  +  c7{r,E)J {f,E)  =  qi{r,E)  +  J dE'  cTsi{r,E'  E)J(r,E') ,  (124) 

where 

qo(r,E)  =  jdQq(r,Q,E)  ,  (125) 

q^(r,E)  =  jdQQq(r,n,E)  ,  (126) 

and,  for  azimuthally  symmetric  scattering 

(Tsoir,E'  -»  E)  =jdajdQ'as(r,E'  E,Q'-Q) ,  (127) 

(Tsi(r,£;'^  J5)  =  JdQjdQ'Q'  Q  as(r,E'  ^  E,Q' ■^)  ,  (128) 

where  =Q'  Q . 

In  the  diffusion  approximation,  qi(r,E)  is  assumed  to  be  isotropic.  Then 
the  current  density,  defined  as 

J(r,i;)sJdQQ^«^(r,Q,i;),  (129) 

can  be  represented  in  terms  of  the  scalar  flux,  ^(r,E) ,  as 
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J{r,E)  =  -D{r,E)V(l>{r,E) 


(130) 


where 


D(f,E)  =  i[CT(f,£)  -7ii(r,E)<75o(r,£)]‘‘ , 


and  the  average  angle  of  scatter,  ni{r,E),  is  found  from 


E 

j^dE"jdQ'jdQQ'-Q<Ts(r,E-^E’',Q'Q) 

= - - - - - 

dE"jdQ'jdQ(7s(r,E-->^E",Q'  Q) 


(131) 


(132) 


Equation  (130)  is  used  instead  (124)  to  eliminate  J(r,E)  from  equation 
(123)  and  gives  the  diffusion  equation 

-W\D{r,E)V<l>(r,E]  +  cT{r,E)<l>{r,E)  =  qQ{r,E) 

/I  qo> 

+^dE'(Jso{r,E' -^  E)<l){r,E'). 

Knowing  the  total  cross  section,  cr(r,£’) ,  the  scattering  cross  section, 
o'so(^>^)  >  and  the  average  scattering  angle,  7/o(F,£’),  we  can  calculate  the 

diffusion  coefficient  and  use  it  in  diffusion  problems.  T-Scat  tracks  all  these 
quantities  (groupwise)  as  is  discussed  below. 

Validity  of  the  MC  Facet  Method  in  the  Diffusion  Limit 

Proper  calculation  of  both  the  scalar  flux  and  current  are  necessary  for 
proper  transport  code  operation  in  the  diffusion  limit.  To  accurately  calculate 
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the  currents,  we  need  the  scattering  cross  sections  to  be  accurate  at  least  to 
the  first  angular  moment  of  scatter.  Further,  the  facet  quadrature  employed 
must  integrate  the  zeroth  and  first  moment  of  the  angular  flux  accurately. 
Conserving  the  first  angular  moment  of  the  scatter  is  described  in  the 


following  section.  Here,  we  perform  a  quick  check  on  the  accuracy  of  the 
Cartesian  quadrature  (with  polecaps). 

The  —  in  equation  (124)  arises  as^^ 

3 

jdQQQ  I 
J  dQ  ~  3 

For  our  quadrature,  we  would  want 

1 

JdQ  JdQ  jdQ  3' 


10  0 
0  10 
0  0  1 


(134) 


To  get  proper  performance  in  the  diffusion  limit,  we  need 

N 

re=l _ _£ 

^  3 

n=l 


(136) 


exactly,  and  similarly  for  rj  and  ^  direction  cosines.  As  long  as 
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JdQQ 


o  - 


JrfQQ 


An,, 


(137) 


this  will  never  be  exact.  Still,  the  values  are  close.  For  the  74  facet 
quadrature  described  earlier, 


N  N 


n=l 


M=1 


N  N 

71=1  71=1 


=  0.332, 


(138) 


and 


^ - =  0.336, 

ZAn» 

71=1 


(139) 


corresponding  to  less  than  a  1%  error  in  the  approximation. 


T-Matrix  Rebalance  to  Conserve  First  Angular  Moment 

Just  as  diffusion  theory  requires  knowledge  of  the  first  angular  moment  of 
the  scattering  cross  section,  for  diffusion-like  problems,  the  T-matrix  should 
accurately  represent  the  angular  structure  of  the  scatter  to  first  order. 
Secondly,  it  would  be  useful  if  the  diffusion  coefficient  could  be  calculated 
directly  from  the  T-matrix,  providing  another  use  for  the  T-matrix  scattering 
libraries. 
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The  facet  tallies  described  earlier  only  recorded  zeroth  order  scatter. 
T-Scat,  however,  also  tracks  the  first  angular  moment  of  the  scatter  into  each 
facet  for  each  energy  group.  Knowledge  of  the  first  angular  moment  within  a 
facet  is  used  to  adjust  the  zeroth  moment  of  the  facet's  neighbors  to  preserve 
the  first  angular  moment  within  the  discrete  ordinate  structure.  This  section 
describes  how  that  is  done. 

For  discrete  ordinates,  a  single  angular  ordinate  represents  a  facet.  If  we 
simply  sum  the  total  number  of  counts  within  a  facet  and  determine  the 
average  (based  on  their  weights),  we  have  a  zeroth  order  approximation  and 
lose  the  information  concerning  the  average  scattering  angle  within  the  facet. 
It  is  possible  within  the  T-Scat  code,  however,  to  track  the  average  scattering 
angle  within  a  given  facet.  As  the  T-matrix  is  calculated,  for  each  count 
within  a  given  scattered  facet  (and  energy),  we  record  the  angular  direction  of 
the  scattered  particle,  Q,- ,  and  its  appropriate  scattering  weight  (equation 

(69)).  The  average  of  the  scattering  direction  (Q )  is  calculated  from 


N 


An'n  _  t=l 
^^g’g  -  II N 


(140) 


and  is  output  from  the  code  along  with  the  zeroth  moment  T-matrix.  This  is 
graphically  represented  in  Figure  15. 
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Figure  15.  Tracking  average  scattering  direction  within  a  facet. 

Process  of  Rebalancing 

There  are  three  general  cases  for  the  polecap  quadrature  that  must  be 
considered; 

1.  Rebalance  among  rectangular  facets  {4-point). 

2.  Rebalance  from  a  rectangular  facet  to  the  polecap  and  a 
neighboring  rectangular  facet  {3-point). 

3.  Rebalance  from  a  polecap  to  two  rectangular  facets  {3-point). 

In  each  case,  two  rules  must  be  followed: 

1.  Conserve  the  average  angle  of  scatter. 

2.  Conserve  the  zeroth  moment  of  the  scatter. 


Case  1;  Rebalance  between  normal  rectansular  facets  (4-Doint). 

For  simplicity  we  discuss  facets  in  the  northern  hemisphere  only. 
Realize  that  the  southern  hemisphere  is  handled  with  a  simple  reflection  of 
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the  theory  presented.  Consider  a  scattered  facet  in  which  the  first  moment  of 
the  scatters  lies  in  the  upper-right  quadrant  of  the  facet,  here  labeled  as  Fi 
(Figure  16). 


Figure  16.  Case  1  rebalance. 

The  azimuthal  angle  is  rotated  so  that  ca  =  0  at  the  centroid  of  Fi.  To 
do  this,  let 

w  =  co-cOj,  (141) 

where  the  index,  j,  corresponds  to  the  facet  being  rebalanced.  With  this 

* 

substitution  wi=u>^=0,  and  W2=W4=  Aco  .  The  polar  variable,  requires 
no  modification.  The  average  direction  of  scatter  is  located  at  (So,(^o) .  In 
terms  of  the  directional  cosines,  the  average  direction  of  scatter  is  (?io>7o>^o)  • 

Let  S  be  the  initial  zeroth  moment  of  the  scatter  found  for  Fi.  After 
rebalancing,  S  is  apportioned  to  the  appropriate  nearest  neighbors  adding  to 
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the  zeroth  moments  of  Fi,  F2,  F3,  and  F4  by  31,82,83,  and  84  respectively. 
The  directional  ordinates  of  the  four  facets  in  terms  of  (S,  or  {jl,  rj,  are 
simply 


n,  =  (o,#a)  =>  =  (142) 

Oj  =  =■  (?'2>’72.^a)  =  [t/1“^A^  cos(A<i)),-(/l-^A^  sin(Affl).^A  j,  (143) 

Qs  =(0,fB)=>  =  (i/i-Ib*  AIb)  (144) 

£24=  (Aib.Ib  )  =>  (?4 >^b)  =  fl/l  “  C08(A®),7i  -  sin(Affl),^B  j.  (145) 


Rule  1:  Conserve  the  average  angle  of  scatter.  To  ensure  that  the  average 
angle  of  scatter  is  conserved  we  set  the  tangents  of  the  original  and 
rebalanced  directional  ordinates  equal  to  each  other  such  that 


Mo  ^  +S2M2+S3M3+S4714 

no  Si^l +S2^2 +53^3 +S4^4  ’ 


(146) 


and 


Mo  ^  ^iMl  +S2M2  +^3^3 
^o  ^I'^l +^2^2 +^3^3 +^4^4 


(147) 


We  now  have  two  of  the  four  equations  necessary  to  solve  the  rebalance 
distribution.  Conservation  of  the  zeroth  moment  with  a  self-imposed 
constraint  (described  below)  provides  the  other  two. 
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Rule  2:  Conserve  the  zeroth  moment  of  the  scatter.  A  4-point  rebalance 


allows  an  extra  degree  of  freedom  in  determining  how  the  zeroth  moment  is 
shared  between  the  four  facets.  Here  we  choose  a  bi-linear  approach.  Of  the 
total  amount  of  S,g  is  a  fraction  redistributed  to  the  right,  and  his  a  fraction 
redistributed  up,  so  that 

si=(l-^)(l-MS,  (148) 

S2=^(l-/i)S,  (149) 

Ss={l-g)hS,  (150) 

and  s^=ghS.  (151) 

Summing  these  equations  demonstrates  zeroth  moment  conservation 

Si  +  §2  "t"  ~  S .  (152) 

Thus  equations  (146)  through  (151)  constitute  six  equations  in  six  unknowns 
(si, 82,53,84,^,  and  h). 

Mathematica  was  used  to  solve  this  system  of  six  equations  in  six 
unknowns.  The  resultant  rebalance  equations  are 

a,(aio  ,^0)  =  S  sin(A«  -  5o) ,  (163) 

=  S  8m(5o)  ,  (154) 

.3(So.«o)  =  -S^%^%^sin(A<»-So),  (155) 

D[Q}Q,gQ) 
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(156) 


where, 


N(Sk  .So  .lo)  =  ^  4Vi-4csc(^^)  -  4(i^l-4k  cot(^^)sec(^'^  -  So) 


,(157) 


where  4k=4A'^4B<  and 

D(So  ,4)  =  4t,(Ji-4A  -  Vl-d)cos(A^)  +  (4  b  - 1,4  )Vl-d  co8(A^  -  So)  .(158) 

Case  2:  Rebalance  from  a  rectansular  facet  involving  the  polecap. 

The  next  case  to  consider  is  the  rectangular  facet  to  polecap  rebalance. 

Again,  the  first  moment  lies  in  the  upper,  right  quadrant  of  a  rectangular 
facet,  directly  above  which  is  the  polecap.  The  geometry  for  this  case  is 
shown  in  Figure  17. 


Figure  17.  Geometry  for  3-point  rebalance. 
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are 


The  directional  ordinates  of  the  three  facets  in  terms  of  (S,  or  {Ji,  rj, 


(159) 


Q2  =  =  ^  cos(Aiy),-Jl-,^  ^sin(A6>),,^  j,  (160) 

and,  Q3  =  (ArA,l)=>(^3,^3,^3)  =  (0,0,l).  (161) 

As  before,  we  set  the  tangents  of  the  directional  ordinates  equal  to  each 
other  such  that 


Mo  si/zi  +  S2//2  +  S3JU3 
Vo  Sl^l  +  ^2772  +  S3  73 


(162) 


and 


Mo  _  ^iMi  +  S2//2  +  S3/73 
Si^l +S2^2 +S3^3 


(163) 


We  now  have  two  of  the  three  equations  necessary  to  solve  the  rebalance 
distribution.  Conservation  of  the  zeroth  moment  provides  the  other: 


S  —  +  S2  ®3  • 


(164) 


The  resultant  rebalance  equations  are 


/-  e  X  oVl“'^o  csc(Aa)) 

Si(6>o.^o)  =  ‘S — ^  ^ — sin(A©-GJo) . 


):  \  cVl-^O  CSC(A®)  . 


(165) 


(166) 


74 


®3(^0>^o)  -  ^ 


^0  ^1-#^  -  Wl-^0  cos(^^  -So)  sec(^^) 


^(Sq.^o) 


(167) 


where, 


i)(So .^o)  =  r - (#  - l)^/wi  cos(A^  3„)sec(A^) 


(168) 


Case  3:  Rebalance  from  a  polecap  to  two  rectansular  facets. 

Rebalancing  from  the  polecap  to  its  neighbors  is  done  in  exactly  the  same 
way  as  the  case  2  rebalance.  The  concept  of  “nearest  neighbor”,  however, 
leads  to  the  conundrum  -  which  neighboring  facets  are  used  in  the  rebalance? 
A  straightforward  approach  is  to  say  that  a  first  moment  slightly  offset  firom 
the  pole  shares  directly  with  its  nearest  neighbors  (along  the  same  meridian), 
and  only  a  3-point  rebalance  is  performed.  Another  argument  is  that  such  a 
slight  offset  should  be  distributed,  in  an  appropriate  manner,  with  all  the 
pole's  neighbors  in  the  first  moment's  hemisphere.  The  first  method  does  not 
take  into  account  that  the  scattering  into  a  facet  may  not  be  uniform. 

Suppose,  for  example,  that  scattering  into  the  pole  is  bi-modal,  such  that  very 
few  neutrons  scattered  directly  into  the  center  of  the  pole,  while  a  large 
number  scattered  to  one  side  and  the  other.  The  resultant  first  moment 
would  be  very  near  the  pole,  and  when  rebalanced,  might  share  with  a  facet 
perpendicular  to  the  direction  the  neutrons  truly  scattered. 
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To  account  for  this  effect  in  a  consistent  manner,  the  polecap  is  itself 
divided  into  the  same  number  of  longitudinal  facets  as  the  rest  of  the  sphere 
as  is  shown  below. 


Figure  18.  Division  of  the  polecap  into  sub-facets. 


The  first  moment  of  the  scatter  is  tracked  within  each  of  these  sub-facets. 
It  would  be  nice  to  use  these  sub-facet  first  moments  to  rebalance  directly 
with  the  nearest  neighbors.  Unfortunately,  you  can't  directly.  Neutrons 
scattered  isotropically  into  the  pole  would  result  in  first  moments  in  each  of 
the  sub-facets  that  were  offset  from  the  pole.  As  such,  each  would  share  with 
its  neighbors,  when  in-fact  no  sharing  should  take  place  at  all. 

Instead,  our  approach  is  to  examine  each  sub-facet  with  its  antithetic  twin 
(the  sub-facet  directly  across  the  pole).  The  two  first  moment  vectors  are 
summed  to  obtain  a  new  first  moment  direction  for  the  pair  , 


76 


Q 


new 


(rpn'n  ^n'n\  ,  (rpn*n 

^g'g  ^^g^gj'^y^g'g 


twin 


(169) 


the  distribution  of  scattered  neutrons  into  the  pole  is  isotropic  the  result 
is  a  first  moment  that  aligns  with  the  pole.  In  some  cases  the  new  first 
moment  may  not  point  within  either  of  the  twin's  facets.  The  facet  the  new 
first  moment  lands  in  is  recorded  for  each  facet  pair.  The  magnitude  of  the 
contribution  given  as 


rnn'n{new)  _  rpn'n  , 
^g'g  -^g'g^ 


(170) 


The  new  superfacet  is  used  in  a  3-point  rebalance  with  its  neighbors  as 
described  above. 


Results  of  the  First  Moment  Rebalance, 

In  this  section  we  look  at  the  meaning  of  the  rebalance  equations  for  the 
various  cases.  This  should  shed  some  light  on  to  what  the  equations  are 
actually  doing.  In  later  chapters,  we  will  see  how  the  rebalance  impacts  the 
T-matrix  in  a  global  sense. 


Expectations  for  cell  rebalance.  On  a  rectangular  grid,  as  shown  in  Figure 

19,  consider  the  average  angle  of  scatter,  Q ,  to  lie  in  the  upper  right 
quadrant  of  a  facet  (Fi),  with  neighbors  F2,  F3,  and  F4  as  shown. 
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Figure  19.  Rebalancing  between  4-points. 

If  Q  were  to  point  in  the  same  direction  as  the  directional  ordinate  of  Fi, 
we  would  expect  Fi  to  receive  one  hundred  percent  of  the  contribution  of  the 

zeroth  moment  of  the  scatter,  and  the  neighbors  to  receive  nothing.  If 

^  ^  _ 

Q  =  Q.g  directly  between  Fi  and  Fs  as  shown  in  Figure  19,  then  these  two 
facets  should  share  the  zeroth  moment  while  the  others  receive  nothing. 
Similarly,  if  Q  =  directly  between  Fi  and  F2  as  shown  in  Figure  19,  then 
these  two  facets  should  share  the  zeroth  moment  while  the  others  receive 
nothing.  Finally,  if  Q  =  ,  we’d  expect  all  four  facets  to  share  equally  in  the 

zeroth  moment. 

Figure  20  shows  what  happens  as  Q  sweeps  from  Q  =  Qj  to  Q  =  . 

Here  Q  lies  on  the  equator  =0)  between  Fi  and  F2  with  no  polar 

component.  For  =  0 ,  Q  =  and  no  rebalancing  takes  place.  As  Q 
moves  toward  >  F2’s  fraction  of  the  zeroth  moment  increases.  The  most  it 
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can  receive  is  half  the  total  contribution  to  the  original  amount  found  in  Fi. 
In  this  case  F3  and  F4  receive  no  contribution. 


Rebalance  Distribution 
Equatorial 


Figure  20.  Rebalance  fractions  between  Fi  &  F2  with  no  polar  component. 

Care  must  be  exercised  when  applying  flat  surface  thinking  to  the  curved 
surface  of  a  sphere.  The  lines  tend  to  distort  from  our  linear  expectation  of 
the  flat  surface  as  one  moves  up  the  pole.  Consider  a  sphere  divided  with  six 
latitudinal  and  twelve  longitudinal  facets  with  two  polecaps  (Figure  6). 

That’s  6x12+2=74  facets  of  equal  solid  angle  on  the  unit  sphere.  Let  Fi’s 
directional  ordinate  lie  on  the  polar  cosine  =  0.243.  The  corresponding 

ordinate  for  F3  lies  on  the  polar  cosine  ^3  =  0.405.  As  before,  let  Q  run  across 
the  middle  of  Fi  from  Q  =  to  Q  =  .  Then  Q  has  a  polar  cosine 

4b  =  0.243  with  Sp  e[0,Afi)/2] .  The  results  are  shown  in  Figure  21. 
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Rebalance  Distribution 
Facets  1  &  2 


Figure  21.  Applying  flat  thinking  to  curved  surfaces:  4-Point  Rebalance. 

One  can  see  that  there  are  negative  fractions  going  to  F3  and  F4,  which  of 
course,  is  a  mathematical  artifact.  We’ve  overstepped  the  domain  of  the 
solution  by  trying  to  follow  the  latitudinal  arc  between  Fi  and  F2  when  we 
can’t  rebalance  past  the  great  circle  arc  on  which  both  and  Q2  lie  (see 
Figure  22). 


Rebalance  Distribution 
Facets  3  &  4 


cdo/Ao)  across  facet 


Figure  22.  Great  circles  define  the  shaded  area  within  which  the  directional 

ordinate  must  lie. 
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The  problem  is  solved  in  the  following  way.  When  determining  the 
“nearest  neighbor”  facets,  one  needs  to  verify  if  Q  is  above  or  below  the  great 
circle  arc  on  which  both  Qj  and  Q2  ’^he  answer  is  found  from 


^  /V 

equation  (158).  If  the  polar  cosine  of  is  and  the  polar  cosine  of  Q 

within  Fi  is  ,  then  the  nearest  neighbor  is  determined  by  comparing 


COS 


Ao) 


with  cos 


Aco  ~ 


where,  as  before. 


coq  e[0,  A<u/2].  If 


cos 


Aco 


cos 


Ao) 


-CO, 


(171) 


then  F3  and  F4  are  the  neighbors  directly  above  Fi  and  F2  (in  the  northern 
hemisphere).  If  not,  then  the  neighbors  lie  below  Fi  and  F2. 

This  same  requirement  applies  for  both  3-point  and  4-point  rebalancing. 
Using  the  same  74  facet  quadrature  described  above,  the  polar  cosine  of  the 

directional  ordinate  of  Fi  (just  below  the  polecap)  is  =  0.905.  Sweeping  Q 
from  Sq  =  0  to  Aco  12  with  ^0=  0.95  across  Fi  the  rebalance  fractions  are 
shown  in  Figure  23. 
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Rebalance  Distribution  Polecap 


Figure  23.  Rebalance  at  the  pole.  (<^i  =  0.905,  <§)  =  0.95) 

As  expected,  for  this  intermediate  polar  cosine,  the  fractions  tend  to  be 

split  50/50  as  Q  moves  across  Fi  toward  F2.  Notice  that  the  polecap  fraction 
tends  to  droop  as  the  sweep  tends  to  Ao  /  2 .  This  is  expected  behavior 
because  of  the  curved  surface  of  the  sphere. 

General  Summary  of  First  Moment  Conservation 

The  process  of  the  rebalancing  code  is  summarized  in  Algorithm  2. 
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Algorithm  2.  First  Angular  Moment  Conservation 


Set  desired  facet  quadrature,  total  facets  =  N. 

Set  desired  group  structure,  total  groups  =  G. 
Calculate  directional  ordinates, 

Dog'=l..G 
Do  ii'=l..N 

BalancedT  =  0  Array  dimension  (N  x  G) 

Read  unbalanced  Tgg  ( g=\..G;  n=\..N) 


Read  facet  average  angular  scatters,  Qg>g  ( g=l..G;  n=l..N) 
For  poles 

Find  new  direction  for  antithetic  twins: 


Determine  new  landing  facet:  n(new). 
Calculate  magnitude  at  this  direction: 


rpn*n{new)  _ 
g'g 


twin 


Save  values  as  separate  contributors  for  later  rebalance. 

For  each  rectangular  facet  and  superfacet  contribution 
Determine  appropriate  nearest  neighbors 
Determine  3-point  or  4-point  rebalance 
Adjust  for  walks  off  left  or  right  side  of  Cartesian  map 
For  appropriate  neighbors  i=l..M  (M  =  3  or  4) 

Determine  appropriate  portion  for  self  and  each  neighbor,  s(i) 

BalancedT^, =  BalancedT^, +  s(i) 


Apply  rotational  symmetry  to  BalancedTgg  : 

For  azimuthally  symmetric  facets  nl  and  n2 

BalancedTgg^  =  {^BalancedT gg^  +  BalancedT gg^^  12 

BalancedTgg^  =  BalancedTgg^ 

Write  BalancedTgg  iox  {g=l..G;  n=l..N) 

EndDo:  g’ 

EndDo:  n 
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The  algorithm  includes  an  adjustment  for  azimuthal  symmetry.  Because 
the  calculation  of  the  facet  contributions  is  performed  with  MC  draws, 
calculation  of  scatter  azimuthally  to  the  left  or  right  of  the  incident  facet  may 
differ  dependent  on  the  number  of  draws  taken.  Once  facet  rebalance  has 
taken  place,  the  azimuthally  symmetric  facets  (now  centered  on  their 
directional  ordinates)  are  simply  summed  with  half  the  total  going  to  each 
facet. 

To  conserve  the  first  angular  moment  of  the  scatter,  we’ve  developed  a 
method  that  tracks  the  average  angle  of  the  scatter  within  the  cell  and  then 
rebalances  the  result  to  the  appropriate  nearest  neighbors.  We  would  expect 
this  to  improve  our  solution  to  transport  problems  that  are  diffusion-like  in 
nature.  While  the  transport  code  we  use,  TETRAN,  is  unable  to  perform 
diffusion  like  problems,  later  we  will  see  the  effects  of  rebalancing  on  the 
T-matrix  itself,  and  then  look  at  how  rebalancing  impacts  the  solution  of  the 
transport  equation  for  fast  neutron  problems. 
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8.  Code  Development,  Capabilities  and  Operation 

The  creation  and  use  of  a  T-matrix  library  is  a  multi-step  process.  It 
begins  with  the  ENDF  data  and  ends  with  results  from  a  transport  code  that 
uses  the  T-matrix  library.  Such  results  are  often  compared  among  codes  for 
validation  or  benchmarking.  In  this  chapter  we’ll  briefly  look  at  the 
validation  process,  discuss  the  libraries  necessary  to  create  both  the  T-Scat 
library  and  the  libraries  used  to  generate  the  spherical  harmonics  cross- 
sections.  We’ll  then  take  a  close  look  at  the  steps  required  to  get  the  data 
generated  from  T-Scat  into  a  form  that  the  transport  code  can  read. 

From  Continuous  to  Discrete:  The  Path  for  Benchmarking. 

To  determine  the  validity,  and  often  accuracy  of  a  discrete  ordinates,  Sj^ , 
code,  programmers  will  typically  compare  results  of  test  problems  with  Monte 
Carlo  transport  techniques  like  MCNP"^-  24.  Often,  the  argument  is  that  MC 
methods  do  not  suffer  from  the  unknown  effects  of  angular,  group  and  spatial 
discretization  found  in  methods. 

There  are,  inevitably,  artifacts  from  the  discretization  schemes  employed 
in  methods,  just  as  there  are  artifacts  resultant  from  variance  reduction 
techniques  employed  by  Monte  Carlo  methods.  Because  the  two  methods  are 
so  very  different  at  their  core,  it  is  unreasonable  to  assume  that  an  Sn 
method  would  converge  exactly  to  a  Monte  Carlo  solution.  The  S]^  solution 
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may  approach  the  MC  solution,  giving  the  user  a  sense  of  well  being,  but  was 
the  convergence  primarily  because  of  a  reduction  in  the  spatial  grid,  angular 
quadrature,  group  structure  or  particular  combination  of  the  three? 

We  propose  a  methodology  to  examine  such  questions,  and  in  so  doing 
introduce  a  new  type  of  transport  code,  developed  at  AFIT,  that  takes 
advantage  of  strictly  positive  T-matrix  data. 

The  Triad  of  Discretization 

One  path  from  a  purely  continuous  method,  like  MCNP,  to  a  discrete 
ordinates  method,  like  TETRAN''  (used  in  this  research),  is  represented  in 
Table  2. 


Table  2.  Discretization  path  from  MCNP  to  Discrete  Ordinates. 


Transport  Code 

Space 

Angle 

Energy 

MCNP 

C 

C 

C 

MCNP-MG 

C 

C 

D 

MCSN 

c 

D 

D 

Discrete 

Ordinates 

D 

D 

D 

C  =  Continuous;  D  =  Discretized 

As  discussed  earlier,  MCNP  is  a  Monte  Carlo  Transport  code  developed  by 
LANL.  It  has  an  option  to  use  standard  Legendre  polynomial  based  multi- 
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group  libraries,  and  calls  this  option  MCNP-MG.  MCSN  is  a  new  transport 
code  developed  at  AFIT^s. 

The  MCSN  method  uses  a  continuous  Monte  Carlo  approach  to  transport 
the  particles  in  space.  The  transport,  however,  is  discretized  in  energy  by 
group,  g,  and  in  angle  by  facet,  n,  represented  by  the  directional  ordinate, 

.  For  each  material,  incident  facet,  n',  and  group,  g',  MCSN  uses  the 
T-matrix  to  create  a  cumulative  distribution  function  of  scattered  facet  and 
group,  F(  n ,  g;  mat,  n' ,  g').  Given  the  material  and  an  incident  facet  and 
group,  MCSN  inverts  F(  n ,  g;  mat,  n',  g').  to  find  the  scattered  facet  and 
group.  It  then  transports  the  scattered  particle  again  in  direction  with  a 
MC  sampling  continuous  in  space. 

This  research  employs  Tetran  as  the  discrete  ordinates  method  of  choice. 
This  discrete  ordinates  code  runs  on  geometries  that  are  spatially  discretized 
with  unstructured  tetrahedral  meshes.  It  can  use  either  standard  Legendre 
polynomial  based,  multigroup  libraries  to  generate  T-matrices  internally 
(potentially  negative)  or  it  can  read  T-matrices  directly  calculated  from 
T-Scat  (guaranteed  non-negative). 

We  emphasize  here  that  the  creation  of  the  non-negative  T-matrix  has 
made  MCSN  possible.  When  sampling  groups  and  ordinates,  the  scattering 
cross  sections  constitute  a  discrete  distribution,  f{n,  g).  It  is  meaningless  to 
have  negative  values  as  would  be  created  from  a  spherical  harmonic 
approximation.  If  they  are  included,  the  cumulative  distribution,  F(n,  g),  is 
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not  monotonic  and  cannot  be  inverted  to  find  g  and  n.  We  can  now  see 
directly  the  effects  of  spatial  discretization  by  comparing  TETRAN  results  to 
MCSN.  We  can  see  directly  the  effects  of  angular  discretization  by  comparing 
results  between  MCSN  and  MCNP-MG.  Finally,  group  discretization  effects 
can  be  examined  by  comparing  MCNP-MG  with  MCNP.  While  the 
exploration  of  discretization  effects  in  various  transport  codes  is  not  the 
thrust  of  this  study,  this  research  has  now  made  such  exploration  possible. 

In  subsequent  chapters,  we  present  some  early  results  comparing  such 
methods.  It  is  certainly  an  area  for  greater  research. 

Cross  Section  Libraries:  Getting  ENDF  data  into  the  transport  code. 

A  number  of  libraries  have  combined  experimental  observations  with 
theoretical  models  to  achieve  a  reasonable  database  of  the  collision  process. 
Many  of  these  libraries  are  available  on  the  web^s.a?  data  prepared  for 
MCNP  comes  from  the  Evaluated  Nuclear  Data  Files  or  ENDF.  The  latest 
release  of  this  library  is  ENDF/B-VI^s.  All  our  calculations  begin  with  this 
cross-section  library.  It  should  be  noted  that  the  accuracy  of  these  cross 
sections  can  vary  considerably.  It  is  important  therefore,  that  for 
comparisons  between  transport  methods,  such  methods  use  or  originate  fi'om, 
whenever  possible,  the  same  cross  section  library. 
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NJOY:  Creating  Continuous  and  Multi-Group  Data. 


There  are  also  a  number  of  codes  that  may  be  used  to  convert  ENDF  data 
to  a  form  used  by  different  transport  methods.  Again,  because  we  will  later 
compare  with  MCNP,  we  use  a  code  called  NJOY,  developed  by  LANL,  to  get 
us  one  step  closer  to  the  libraries  we  need.  NJOY  contains  a  number  of 
subprograms  that  read  in  the  ENDF  data  and  convert  it  to  either  continuous 
energy  cross  section  libraries  (as  used  by  MCNP)  or  to  create  the  Legendre 
coefficients  used  to  calculate  the  multi-group  spherical  harmonics.  The 
continuous  energy  format  is  known  as  an  ACE  file.  The  ACE  format  is  read 
by  MCNP  to  perform  transport  and  also  by  T-Scat  to  create  the  T-matrices. 
The  output  form  of  the  multi-group  data  is  known  as  a  MATXS  library. 
Another  code  developed  by  LANL,  TRANSX^^,  reads  the  MATXS  library  and 


Figure  24.  Flow  of  the  Cross-Section  Data  to  the  Transport  Codes 
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creates  a  new  library  (CARD)  that  can  be  used  by  standard  transport 
codes  (using  spherical  harmonics  to  approximate  the  scattering  source). 
Figure  24  depicts  the  flow  of  cross-section  data  through  the  various  processes. 
The  codes  that  perform  the  processing  are  shown  in  the  rounded  boxes.  The 
darkened  boxes  represent  those  codes  and  libraries  created  at  AFIT. 


Group  Structures 

The  NJOY  processing  system  allows  the  user  to  select  a  number  of 
predetermined  group  structures.  In  this  research  we  use  two,  LANL-30  and 
VITAMIN-J.  Both  structures  can  be  found  in  the  appendix.  The  LANL-30 
ranges  from  1.39E-04  eV  to  17.0  MeV  and  contains  30  groups.  The  group  size 
is  defined  using  the  quantity  called  lethargy,  u,  where 

u(£;)  =  ln^.  (172) 


When  using  lethargy  to  describe  the  size  of  the  group,  we  define  the  group 
lethargy  as 


Am.  sln^-ln-^ 

^  Eg  Eg-i 


=  ln- 


E 


g-i 


E 


g 


(173) 


A  group  size  of  one  lethargy  implies  that 


Am^ 


=  ln- 


E 


g-\ 


E 


=  1, 


g 


or 


(174) 
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(175) 


For  the  LANL-30  group  structure  (Figure  25  ),  the  group  size  at  low 
energy  is  one  lethargy  and  at  higher  energies  %  and  ^4  lethargy  groups  are 
used.  It  is  essentially  a  survey  group  structure,  allowing  the  researcher  to 
get  a  rough  idea  of  the  problem  at  hand.  This  group  was  selected  because  it 
is  the  group  structure  of  MATXS-10,  a  standard  MATXS  library  used  to 
develop  the  Legendre  coefficients  for  transport  comparisons. 

The  Vitamin-J  structure  is  an  evolution  of  an  older  library  used  for 
shielding  applications.  The  energy  bounds  range  from  l.OE-05  eV  to 
19.64  MeV  with  a  total  of  175  energy  groups. 


Lethargy  by  Group  Structure 


Figure  25  Group  size  in  lethargy  for  selected  group  structures 
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This  group  was  selected  to  help  test  the  anisotropy  that  a  finer  group 
structure  can  create.  Further,  Tetran  was  developed  to  employ  the  EC 
method,  ideal  for  shielding  applications.  Taking  advantage  of  this  synergy, 
some  deep  penetration  problems  are  examined  as  part  of  this  research.  Of 
course,  any  group  structure  can  be  used  with  any  of  the  codes,  including 
T-Scat. 

T-Scat’s  Current  Capabilities 

Here  we  outline  the  steps  required  to  get  from  the  ACE  format  to  a  library 
that  is  used  by  the  transport  code  TETRAN.  There  are  three  basic  steps, 
each  representing  separate  codes  developed  for  this  research. 

1.  T-Scat;  Performs  MC  integration  to  calculate  the  T-matrix. 
Calculates  total  cross-section  and  determines  first  angular 
moments  for  each  scattered  facet  and  group. 

2.  T-Balancer:  Rebalances  the  uncompressed  T-matrix  to  conserve  the 
first  angular  moments  of  the  scatter. 

3.  T-Mixer:  Decompresses  and  combines  the  libraries  for  the  problem 
at  hand  (group  material  and  mix  specific).  Converts  microscopic 
data  to  macroscopic  data. 

Calculating  the  T-matrix  with  T-Scat. 

As  discussed  earlier,  T-Scat  evaluates  the  T-matrix 
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by  reading  microscopic  cross-section  data  directly  from  an  ACE  file  and 
computing  the  T-matrix  for  the  given  group  and  facet  structure. 

Scattering  Laws.  The  currently  supported  scattering  laws  (following  the 
MCNP  formalism)  are: 

1.  Elastic. 

2.  Discrete  level  inelastic. 

3.  Continuum:  evaporation  spectrum. 

It  would  be  straightforward  to  incorporate  the  other  scattering  laws 
including  multiplying  (n,2n)  and  fission  reactions.  T-Scat  is  both  downscatter 
and  upscatter  capable.  For  a  complete  list  of  such  laws  see  Appendix  B. 
Symmetries  Applied.  Because  of  the  quadrature  used,  one  can  take 
advantage  of  the  available  symmetry  in  the  scattering  and  reduce  the 
number  of  calculations  that  need  to  be  performed^^.  For  each  incident  group, 
only  the  incident  facets  in  the  northern  hemisphere  along  the  same  meridian 
need  be  calculated.  For  each  incident  facet  we  tally  to  the  appropriate 
scattered  groups  and  facets  on  the  entire  sphere.  The  other  incident  facet 
calculations  are  either  reflections,  rotations,  or  both  of  this  same  calculation. 
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Figure  26.  74'l’acet  Structure.  All  incident  facets  are  rotations  or 
reflections  of  the  shaded  regions. 


Figure  26  shows  the  incident  facets  for  which  T-Scat  calculates.  The 
arrow  illustrates  a  T(56^23)  transfer  for  a  given  g'^g  (counting  facets  left  to 
right  from  bottom  to  top).  This  same  scatter  is  shown  on  the  unit  sphere  in 
Figure  27.  Such  a  scatter  is  the  same  as  T(21^59)  by  reflection,  !r(52->’19) 
by  rotation,  and  r(16^55)  by  reflection  and  rotation  (as  shown).  The  data  is 
stored  in  this  compressed  form.  It  is  the  job  of  another  code,  T-Mixer,  to 
decompress  the  data  and  prepare  it  for  transport. 

Other  Outputs.  Also  output  during  this  step  are  the  total  group  cross- 
section  and  the  first  moment  of  scatter  for  each  facet.  The  total  cross  section 
is  integrated  directly  from  the  continuous  energy  ACE  file  using  the  same 
energy  weight  structure  as  the  T-matrix.  The  Monte  Carlo  integration  of  the 
total  cross  section  takes  the  form 
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Because  the  average  scatter  data  is  a  set  of  three  values 
(^11  =  (/^n  ’Vny^n ))  ®^ch  Tgg  ,  the  data  file  is  somewhat  larger  than  the 

T-matrix  file  for  Tg-g  alone,  but  once  the  rebalance  is  performed,  the  average 
scatter  data  is  no  longer  required. 


Figure  27.  T'(56-»23)  facet  transfer  on  the  sphere. 

Rebalance  to  Conserve  the  First  Angular  Moment 

The  next  step  is  to  rebalance  the  T-matrix.  For  each  ,  the  rebalance 

program  determines  the  appropriate  nearest  neighbors  and  distributes  the 
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zeroth  moment  appropriately  conserving  the  first  moment  input.  The  output 
is  still  in  compressed  form  (only  for  the  specific  incident  facets  calculated). 
T-Balancer  performs  the  following  tasks: 

1.  Given  the  group  and  quadrature  structure  determine  facet  directional 
ordinates. 

2.  Verify  that  all  input  first  moments  are  found  within  the  appropriate 
facet.  This  is  done  as  an  error  check  on  the  input. 

3.  Determine  facet  nearest  neighbors  (as  discussed  earlier). 

4.  Apply  4-point  or  3-point  rebalance  as  appropriate. 

5.  Assign  rotational  symmetry.  A  facet  transfer  azimuthally  to  the  right 
should  be  identical  to  the  transfer  to  the  left.  Statistically  these 
numbers  converge  with  increasing  draws. 

6.  Write  the  balanced  T-matrix. 

The  rebalance  code  run  time  is  orders  of  magnitude  shorter  than  the  T-Scat 
code  run  time  to  calculate  the  T-matrix.  The  entire  program  could  be 
implemented  as  a  sub-routine  of  the  main  program,  T-Scat. 

Prepare  Library  for  the  Transport  Code  Tetran 

T-Mixer  prepares  the  T-matrix  data  for  use  by  Tetran.  The  way  Tetran 
reads  in  the  T-matrix  dictates  the  data  structure  produced  by  the  mixer.  It 
should  be  noted  that  other  transport  codes  may  have  other  requirements. 
Hence,  the  T-matrix  libraries  are  only  dependent  on  the  angular  quadrature 
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and  group  structure  used.  It  is  the  mixer’s  job  to  prepare  the  data  for  the 
problem  specific  transport  run. 

Tetran  reads  in  the  T-matrix  data  in  the  following  way.  First,  is 

read  for  all  facets  and  all  materials  in  the  problem.  The  code  is  run  until  the 
group  1  angular  flux,  is  found  for  all  incident  angles.  Once  known,  this 

flux  is  used  to  calculate  its  contribution  to  the  scattered  source  for  all  down 
scattered  energy  groups  by  reading  T^g  for  all  materials  and  groups,  ^>1. 

TETRAN  then  proceeds  to  the  next  group  and  performs  the  same  calculation 
sequence,  summing  up  the  contributions  from  all  previous  calculations  to  get 
the  total  scattered  source  for  the  present  group.  The  code  sweeps  down  the 
group  structure  in  this  way,  until  the  final  group  is  calculated. 

TETRAN  does  not  require  the  T-matrix  to  be  positive.  It  is  the  spatial 
quadrature  techniques  used  by  TETRAN  that  dictate  requirements  for  source 
positivity.  TETRAN  employs  two  types  of  spatial  quadrature-linear 
characteristic  (LC)^®,  and  exponential  characteristic  (EC)'^.  The  LC  method 
does  not  require  that  the  calculated  sources  be  positive.  The  EC  method, 
however,  does  require  non-negative  sources.  Good  for  deep  penetration 
problems,  the  EC  method  will  fail  if  a  negative  source  is  calculated. 

The  T-matrix,  built  by  spherical  harmonics,  is  given  as 

(178) 

/=0  m--l 
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which,  in  TETRAN,  is  calculated  as  required  during  the  code  run.  The 
chapters  that  follow  will  give  a  direct  comparison  between  the  T-matrix  as 
built  by  T-Scat  with  those  built  from  spherical  harmonics.  We  will  also 
present  comparisons  of  transport  code  runs  using  both  the  LC  and  EC  spatial 
quadratures. 

The  next  step  is  performed  by  the  code  T-Mixer.  T-Mixer  prepares  the 
T-matrix  for  transport  use.  If  the  transport  problem  looks  at  multiple 
materials  (e.g.,  water-iron-water)  then  the  T-Scat  matrix  for  water  has  to  be 
shuffled  with  the  T-Scat  matrix  for  iron  as  required  by  the  transport  problem. 
Further,  as  was  described  by  the  section  “Select  the  Collision  Nuclide”  on 
page  49,  the  T-matrix  for  water  (H2O)  may  have  already  been  generated  by 
T-Scat  directly  or  must  be  mixed  from  the  separate  matrices  of  hydrogen  and 
oxygen.  Finally  the  ordinate  structure  of  the  T-matrix  must  be  decompressed 
to  allow  for  direct  reading  of  the  ordinate  to  ordinate  transfers  by  the 
transport  code 

The  mixing  code,  T-Mixer,  first  combines  T-matrix  data  for  mixed 
materials  such  as  H2O  (as  required).  When  there  are  M  different  nuclides 

forming  the  material  then  the  microscopic  cross  sections,  Tj ,  are  combined  as 

M 

Ttot  =  T^ifi,  (179) 

i=l 

where  /j  is  the  atom  fraction  of  nuclide  i.  For  H2O  this  is 
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(180) 


'I'h^o  =  +^o  • 


The  T  -matrices  are  then  converted  from  microscopic  data  (barns)  to 
macroscopic  (cm-i)  data  based  on  the  density  and  isotopes  of  the  material.  If 

T^Q^  is  the  microscopic  T-matrix  calculated  from  equation  (179),  then  the  total 
macroscopic  T-matrix,  '^tot  is  found  as 


M 

N 

Tu,t  =  =  Nav/’# - ,  (181) 

ZAifi  Y.Aifi 

•  1  • 


where,  Nav  is  Avagadro’s  number  (0.60221367  cm^/barn-mol),  p  is  the 
density  (g/cm^)  of  the  material,  and  Aj  the  atomic  weight  (g/mol)  of  nuclide  i 
in  the  material.  For  water,  this  is 


M 

'LTifi 

Tnfi  =  IVAvPir - =  (0.6022)(1.0) 

i=\ 


2Th+To 


2(0.999170) +  15.85316’ 


(182) 


The  T-Mixer  then  selects,  reflects  and  rotates  the  T-matrix,  as 
appropriate,  for  each  incident  facet.  The  data  is  then  written  to  file  to  be 
read  in  by  TETRAN.  The  transport  code  developed  at  AFIT,  MCSN,  reads 
the  data  in  the  same  way,  so  no  separate  calculation  is  required  for  that  code. 

We  are  now  in  a  position  to  evaluate  the  new  non-negative  scattering 
matrix  on  two  fronts.  First,  how  does  the  T-matrix,  generated  with  the  MC 
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facet  method  compare  to  what  we  expect  from  the  physics  and  to  its  well- 
estabhshed  spherical  harmonic  counterpart?  Second,  how  do  transport  codes 
perform  using  this  new  T-matrix  particularly  in  the  presence  of  anisotropic 
scattering?  The  following  chapters  answer  these  questions. 
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9.  Calculations:  Investigations  of  the  T-Matrix 


There  are  two  primary  results  of  this  effort.  The  first  is  the  creation  of  the 
T-matrix  for  various  materials  and  compounds  using  the  Monte  Carlo  facet 
method  and  subsequent  rebalance.  The  second  is  the  performance 
improvement  of  transport  codes  using  this  T-matrix  data.  This  chapter 
examines  the  first  of  these  results. 

In  this  chapter,  well  develop  an  understanding  of  how  to  read  the 
T-matrix  data.  As  was  the  case  when  we  examined  in  one  dimension, 

well  see  that  the  structure  of  the  scattering  matrix  affords  insight  to  the 
physics  of  the  scatter.  Well  show  that  the  Monte  Carlo  facet  method  creates 
a  non-negative  T-matrix  that  directly  represents  the  physics  of  the  scatter. 
Next,  the  convergence  and  run  time  of  the  T-Scat  integrator  is  presented. 
Following  this  is  a  comparison  between  strictly  positive  T-matrix  data 
created  by  T-Scat,  and  data  created  by  the  spherical  harmonic  method.  This 
will  demonstrate  the  superior  performance  of  the  MC  facet  method, 
particularly  for  anisotropic  scattering.  The  next  section  investigates  the 
specific  angular  support  failings  of  ordinate-to-ordinate  methods,  and 
discusses  how  the  MC  facet  method  has  no  such  failings.  Finally,  well 
demonstrate  the  effects  of  angular  rebalancing  on  the  T-matrix  data. 
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Visualizing  the  Scattering  Physics— An  Introduction  to  the  Graphics 


This  section  gives  a  quick  sketch  of  the  various  t5q)es  of  graphics  used  in 
this  research  to  analyze  the  T-matrix  data.  The  goal  is  to  see  what  the 
physics  is  telling  us  and  verify  that  what  we’re  seeing  makes  sense. 

As  discussed  earlier,  a  monoenergetic  particle  travels  with  energy  E',  in 
direction  Q'  and  strikes  a  nuclide.  After  the  collision  the  particle  has  a 
scattered  energy,  E,  and  direction  Q .  Here  we  assume  azimuthal  symmetry 
in  the  scatter  so  that  the  new  direction,  Q ,  is  found  in  some  cone  around  the 

yv 

incident  direction,  Q'  (see  Figure  28) 

n 


Figure  28.  Scattering  Cone  After  Collision. 

Placing  this  geometry  on  the  sphere,  the  cone  is  represented  by  a  circle  on 
the  surface  of  the  sphere  indicating  the  possible  new  direction  of  the 
scattered  particle  (Figure  29).  Folding  the  sphere  out  like  a  map  of  the  earth, 
the  circle,  as  expected,  becomes  distorted. 
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Figure  29.  Scattering  circle  of  the  sphere  on  a  Cartesian  map. 

For  elastic  scatter,  the  greater  the  loss  in  energy  the  greater  the 
scattering  angle,  where  =  cos(^2^) .  Figure  30  shows  the  curves  on 

the  Cartesian  map  corresponding  to  a  progressive  increase  in  0^  and  greater 

energy  loss.  The  scattering  to  the  lower  energies  becomes  almost  square-like 
in  outhne. 


Figure  30.  Progressive  increases  in  scattering  angle  on  the  Cartesian  grid. 
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This  is  important  to  recognize  when  examining  scatter  to  lower  groups  on 
Cartesian  maps,  as  discussed  below. 

Of  course,  the  scatter  is  not  at  a  set  angle,  but  has  a  probability 
distribution  f(jjc,  E').  So  the  line  above  is  really  the  peak  of  a  distribution 
that  falls  off  to  either  side.  Further,  we’re  not  considering  particles  traveling 
within  a  single  direction,  Q' ,  but  rather  isotropically  across  an  incident  facet, 
.  Finally,  the  real  transfer  being  considered  is  not  mono-energetic  but 


multigroup.  This,  en  masse,  tends  to  spread  the  curves  above.  By  accounting 
for  the  probability  distribution  with  a  third  axis  perpendicular  to  the 
Cartesian  map,  a  terrain  is  created  that  illustrates  the  probability  of  scatter. 
If  this  is  done  on  a  group-to-group  basis,  the  data  of  the  third  axis  is  found 
directly  from  the  T-matrix  as  T”'” . 


As  an  example,  we  examine  the  scattering  cross  section  of  hydrogen. 


Hydrogen  Differential  Scattering  Cross  Section 
LANL-30  Group  Structure 


Figure  31.  Group-to-group  differential  scattering  cross  section  of 

hydrogen  as  a  function  of  the  lab  frame  scattering  angle. 
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In  Figure  31,  the  differential  scattering  cross  section  as  a  function  of  lab 


angle  is  shown.  These  were  generated  with  T-Scat  where  the  incident 
neutron  is  traveling  in  group  1  of  the  LANL-SO  group  structure.  The  sphere 
was  divided  into  500  latitudinal  facets.  Only  a  handful  of  the  facets  are  ever 
reached  in  the  first  five  groups  (note  scale). 

Looking  at  the  within  group  scatter  15-17  MeV),  it  is  highly  forward 

peaked.  That  is  for  within  group  scatter  at  high  energy,  neutrons  tend  to 
continue  in  their  original  direction  of  travel.  Placing  this  distribution  on  a 
scatter  plot  like  the  Cartesian  map  of  Figure  29  and  plotting  in  three 
dimensions,  we  would  expect  a  cone  to  appear  peaked  at  the  center  of  the  grid 
and  falling  off  to  the  sides.  The  same  is  true  of  a  facet  transfer.  In  the  facet 
quadrature,  particles  travehng  isotropically  within  the  facet  tend  to  stay 
within  the  facet  after  a  (1— >^1)  collision.  The  scattered  distribution  is  also  a 
cone  in  the  center  of  the  Cartesian  grid. 
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Tin'n 
1^1 


Hydrogen  Within  Group  Scatter 


Figure  32.  Three  dimensional  look  at  within  group  scatter  of  on  the 
Cartesian  grid  for  a  74  facet  quadrature. 

Figure  32  was  created  from  T-matrix  data  where  the  T  values  are 
centered  within  each  facet,  at  .  The  T-matrix  facet  structure  is  the 
74  facet  Cartesian  grid  described  earlier.  The  graphing  routine, 

TableCurve  3D3i,  performed  a  bi-quadratic  order  2/2  spline  interpolation  of 
the  T-matrix  data.  While  the  interpolation  may  swing  slightly  negative,  the 
data  is  strictly  non-negative.  The  raw  facet  values  can  be  seen  by  looking  at 
a  bar  chart  of  the  data  given  above  (Figure  33). 
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Figure  34.  Hydrogen  within  group  scatter  showing  the  explicit  facet 
values  of  the  T-matrix  (21  x  40:  762  facets). 


It  is  also  possible  to  look  at  the  data  explicitly  in  three  dimensions  using 
spherical  coordinates  [r,  co,  cos  i(^)  ].  The  angular  coordinates  are  the  facet 
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ordinates  and  the  T-matrix  value  provides  the  r  value.  The  result  is  a 


sharply  forward  peaked  pattern.  The  most  probable  direction  of  scatter  is 
forward  peaked  with  no  chance  of  back  scatter.  Of  course,  for  a  stationary 
nuclide,  elastic  back  scatter  of  hydrogen  in  the  lab  frame  is  not  possible. 

i  I 


Figure  35.  3-Dimensonsional  representation  of  within  group  scatter  for 
hydrogen  using  the  first  group  of  the  LANL-30  structure). 

The  Cartesian  grid  quadrature  allows  the  data  above  to  be  plotted  directly 
using  Mathematica.  The  directional  ordinates  of  the  standard  level 
symmetric  quadrature  (used  with  spherical  harmonic  expansion),  however, 
makes  such  plotting  more  difficult.  To  facilitate  comparison  between  the  two 
methods,  we  draw  the  data  points  directly  with  rays  extending  from  the 
origin.  In  this  way,  the  individual  data  values  of  Figure  35  are  now  plotted 
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three  dimensionally  in  Figure  36.  When  using  illustrations  like  Figure  36  to 
present  spherical  harmonic  data  negative  values  are  set  to  zero. 


Figure  36.  Exphcit  representation  of  T-matrix  data  as  scatter  points 

The  graphics  shown  above  will  be  used  extensively  in  the  rest  of  this 
dissertation.  These  plots  will  be  used  to  help  illustrate  a  number  of  the 
advantages  of  the  MC  facet  method  T-matrix.  In  the  sections  below,  we  will 
answer  the  following  questions: 

•  Is  the  physics  of  the  scatter  being  maintained? 

•  Does  the  code  converge  as  expected  in  a  reasonable  amount  of  time? 

•  How  does  the  MC  facet  method  compare  to  spherical  harmonics? 


109 


•  Is  within  group  scatter  and  next  group  scatter  guaranteed? 

•  Are  groups  skipped? 

•  What  are  the  effects  of  angular  rebalance? 

In  answering  these  questions  we’ll  find  the  MC  facet  method  superior  to 
its  counterparts  in  a  number  of  important  ways. 

Visualizing  the  physics  of  group-to-group  scatter. 

Figure  31  shows  the  group-to-group  scatter  of  hydrogen  in  1-D.  Such  plots 
are  common  in  the  literature®’^^^  jn  this  section,  we  expand  the  standard 
analysis  to  examine  the  group-to-group  T-matrix  transfers  created  by  the  MC 
facet  method. 

As  shown  in  Figure  30,  the  group  structure  directly  determines  the 
allowed  angle  of  scatter.  Combining  this  type  of  contour  plot  with  the  group- 
to-group  scattering  distribution  of  Figure  31,  the  MC  facet  method 
demonstrates  the  same  behavior.  Such  a  contour  plot  is  shown  in  Figure  37. 
As  the  scattered  neutrons  lose  more  and  more  energy  they  scatter  further 
and  further  away  from  the  incident  direction. 
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Figure  37.  Cartesian  contour  of  the  first  five  scatter  groups  of  hydrogen. 

The  plots  we’ve  shown  to  this  point  have  been  for  hydrogen.  Hydrogen  is 
perhaps  the  easiest  nuclide  to  model  because  its  scattering  is  strictly  elastic 
(for  the  collisions  of  interest).  It  is  perhaps  the  most  difficult  nuclide  to  model 
correctly  because  of  its  sharply  anisotropic  behavior.  In  the  center-of-mass 
(CM)  frame,  hydrogen’s  scattering  distribution  as  a  function  of  incident 
energy,  /(jE',//^,^)  .  is  essentially  isotropic.  Hydrogen’s  anisotropy  comes 


from  converting  the  CM  scatter, //cm  >  lab  frame  scatter,  For  heavier 

nuclei,  isotropic  scattering  in  the  CM  frame  corresponds  to  isotropic  scatter  in 
the  lab  frame.  Anisotropy  in  the  lab  frame,  for  heavier  nuclei  results  from 
anisotropic  scattering  distribution  functions,  /(jEJ',//c,„)  ,  which  tend  to  be 


forward  peaked  at  high  energies.  Inelastic  scattering  levels  and  the  group 
structure  used  introduce  other  anisotropic  behavior. 

As  an  example  of  heavier  nuclei,  consider  oxygen  (^^0).  This  nuclide  has 
38  different  inelastic  levels  that  must  be  included  in  the  physics  model.  Each 
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level  has  a  corresponding  Q-value,  for  ranging  between  6.05  and 
18.55  MeV.  Because  of  the  kinematics  involved,  neutrons  striking  oxygen  at 
these  energies  may  have  two  different  scattered  energies  for  a  given  scatter 
angle,  .  The  scattered  energies  are 

+  (183) 

(A  +  1)  I  *•  J 
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and  Qi  is  the  excitation  energy  for  the  considered  inelastic  level.  A  similar 
relation  for  the  incident  energy  is 
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To  illustrate  these  scattering  kinematics,  the  relation  between  incident  and 
scattered  energy  at  given  values  is  plotted  in  Figure  38.  Clearly*,  the  finer 

the  group  structure  used,  the  more  confined  the  allowed  range  of  jXj^  and  the 
more  anisotropic  the  behavior  of  the  group-to-group  scatter. 
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Elastic  Scatter  Regime 
Group  l->2  (LANL-30) 


Figure  38.  Kinematics  for  the  group  l->2  scatter  of 

The  interplay  between  the  scattering  probability  distribution  function, 
f{E',jUc) ,  the  atomic  weight  of  the  target  nuclide,  and  the  group  structure 

used,  must  be  considered  when  examining  anisotropic  behavior.  Figure  39 
shows  the  group -to-group  T-matrices  for  the  scatter  from  the  first  group 
(15-17  MeV)  to  the  first  four  groups  of  the  LANL-30  group  structure. 


113 


T(l-^l) 


T(l-^2) 


T(1^4) 

0.00005 


“0.00005  0  0.00005 


Figure  39  The  scattering  of  Group  one  to  the  first  four  groups  of  the 

LANL-30  group  structure. 

Note  that  the  four  illustrations  have  different  scales.  Relative  to  the 
group  1->1  transfer,  the  l->2  and  l->3  and  l->4  are  progressively  smaller. 
This  is  consistent  with  the  physics  of  the  scattering  process. 

The  collision  of  a  neutron  with  oxygen  is  far  less  efficient  in  slowing  down 
the  neutron  than  hydrogen.  Of  course,  it  is  the  hydrogen  in  H2O  that  makes 
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water  such  a  good  moderator.  Hydrogen  has  a  fairly  uniform  scattering 
distribution  in  the  first  five  groups  (Figure  31).  A  plot  of  the  T-matrix  for 
oxygen  (using  the  same  group  structure  and  a  21x40-»762  facet  quadrature) 
shows  that  the  group-to-group  probability  distribution  falls  off  rapidly  as  a 
function  of  energy  lost. 


Group-to-group  Scatter 


Figure  40.  Group-to-group  T-matrix  of  on  the  Cartesian  map 

(IjANL-30  Group). 

While  the  (l->2)  scatter  of  Figure  40  seems  smooth,  examination  of  the 
scatter  in  Figure  39  shows  that  it  is  highly  anisotropic.  A  close-up  of  the 
scatter  on  the  Cartesian  plane  demonstrates  that  there  is  a  significant 
amount  of  structure  to  the  scatter,  which  can  not  be  matched  with  a  strictly 
positive  spherical  harmonic  expansion. 
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Oxygen-16  (1'->2) 


Figure  41  Detail  of  (l->2)  T-matrix  scatter  of  on  the  Cartesian  map. 

All  the  scattering  data  presented  so  far  was  created  with  T-Scat.  The  data 
is  strictly  non-negative  and  tends  to  behave  in  a  manner  consistent  with 
what  is  known  of  the  scattering  physics.  Taken  together,  the  data  provides 
some  confidence  that  what  T-Scat  is  doing  is  appropriate.  The  next  step  is  to 
see  how  well  the  facet  integrator,  T-Scat,  is  working.  Particularly  we  want  to 
ensure  such  error  tracks  with  standard  MC  integration.  It  is  also  important 
to  know  if  the  integration  is  practicable.  We’ll  show  that  a  reasonable  error 
can  be  achieved  in  a  reasonable  amount  of  time,  and  provide  some  ideas  for 
reducing  such  error  even  further.  Finally,  a  comparison  of  the  results  of  the 
MC  facet  method  with  that  of  the  spherical  harmonics  method  is  presented. 
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Convergence,  Error  Tracking  and  Calculation  Time 

T-Scat  outputs  a  number  of  files  that  can  be  used  to  study  the  processes  at 
hand.  First,  is  a  data  file  that  describes  the  group  and  quadrature  structure 
used.  It  includes  the  materials  being  investigated,  the  scattering  laws 

applied  and  tracks  the  convergence  for  the  transfer.  For  each  group,  the 
total  cross  section,  cr^,  is  calculated  and  output.  For  each  facet  the 
directional  ordinate,  ,  is  output.  For  each  g'->g  and  n'->n  transfer  T-Scat 
outputs  the  integral  value:  Tg,g  ,  the  angular  first  moment,  ,  the  counts 
in  each  facet,  M”-” ,  and  the  error  for  each  facet,  . 

In  performing  the  integration,  T-Scat  doubles  the  number  of  scattering 
draws  for  each  iteration.  The  code  terminates  when  the  maximum  number  of 
draws  is  reached,  or  the  error  for  one  of  the  facet  values  is  less  than  the 
selected  tolerance. 

After  M  draws,  the  error  for  each  facet  is  calculated  as 
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where  ^  was  defined  in  equation  (70),  and  corresponding  to  a  95% 

confidence  that  the  value  for  that  transfer  lies  within  ^of  the  estimated 
valuers. 
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Figure  42  shows  the  value  of  the  1-^1  group  and  facet  transfer  for  oxygen 
dependent  on  the  number  of  draws  (LANL-30  group,  21  x  40  facets). 


Convergence  of  the  Oxygen  T-Matrix 


Number  of  Draws 

Figure  42.  T  value  of  the  oxygen  1->1  group  and  facet  transfer. 

Expanding  the  ordinate  axis  and  placing  the  error  bars  on  the  data  using 
equation  (186)  above,  Figure  43  shows  that  as  the  number  of  draws  increases 
the  newly  calculated  T  value  is  within  the  95%  confidence  level  of  the 
previous  estimate.  Figure  44  plots  the  fractional  error  and  shows  that  it 
agrees  well  with  the  expected  1/ -Jn  convergence. 
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Figure  43.  T  value  of  the  oxygen  transfer  with  error  bars. 


Figure  44.  Fractional  Error  matches  convergence. 
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Oxygen  (l->2)  640k  Draws 


Oxygen  (l->2)  1.28M  Draws 


Oxygen  (l->2)  10.2M  Draws 


Oxygen  (l->2)  41. OM  Draws 


Figure  45.  Increasing  the  number  of  draws  on  the  ^^0(1^2)  transfer  using  a 

21x40->762  facet  quadrature. 

While  20  million  draws  provide  better  than  half  a  percent  error  in  the 
most  likely  transfer  for  oxygen,  this  is  certainly  not  the  case  for  the  least 
probable  scatters  in  the  transfer  matrix.  Recall  that  the  peak  of  the  i®0(l->2) 
transfer  was  about  50  times  smaller  than  i®0(l^l)  transfer.  A  look  at  the 
Cartesian  map  of  the  T-matrix  for  this  scatter  reveals  the  effect  of  increasing 
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the  number  of  draws  (Figure  45).  For  a  small  number  of  draws  the  shape  of 
the  scatter  transfer  is  discernable  but  not  distinct.  At  aboutlO  miUion  draws 
the  essential  shape  of  the  scatter  is  preserved.  Of  course  this  is  for  a 
21  X  40  ->  722  facet  quadrature  coupled  with  the  LANL-30  group  structure. 
Different  quadratures  will  produce  different  results  dependent  on  the  group 
structure  used. 

Still,  by  its  nature  T-Scat  samples  the  T-matrix  in  those  most  probable 
scatter  regions.  More  particles  land  in  the  areas  they’re  supposed  to.  Thus 
the  code  performs  well  for  anisotropic  scattering  where  the  scattered 
neutrons  are  confined  to  a  small  portion  of  solid  angle.  The  error  changes 
across  the  T-matrix.  Figure  46  shows  the  160(1->1)  transfer  while  Figure  47 
shows  the  absolute  error  of  that  transfer  for  5  million  and  then  20  million 
draws.  The  absolute  error  is  greatest  where  the  integral  is  greatest.  Of 
course,  by  sampling  more  the  error  is  reduced.  While  the  absolute  error 
appears  worse  in  the  peak  regions,  the  fractional  error  is  much  less.  Figure 
48  shows  a  log  plot  of  the  error  for  the  within-group  scatter. 
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O  (1->1)  Transfer 
40  million  draws 


Fractional  Error  ^®0{  1->1) 
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Figure  48.  Fractional  error  of  i®0(l->l)  transfer  (20  million  draws). 

A  wholly  different  look  at  error  is  to  examine  what  happens  when  the 
scattering  is  forced  to  be  one  group,  isotropic.  In  this  case  all  the  facets 
should  have  the  same  value.  For  a  10  x  16^130  facet  structure,  we  plot  the 
fractional  error,  Sf,  where 


(187) 


as  a  function  of  the  number  of  draws  in  Figure  49.  The  fractional  also  tracks 
with  the  expected  1/ Viv  convergence. 
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Fractional  Error  v.s.  Draws 


Number  of  Draws 


Figure  49.  Fractional  error  of  one  group,  160  facet  isotropic  scatter. 

One  variance  reduction  technique  conceived  but  not  implemented  in  the 
code  is  the  splitting  of  the  scattered  azimuthal  (coi)  draw.  Because  the 

scatter  is  azimuthally  symmetric  about  the  scattered  angle,  a  split  of  100  or 

360  neutrons  could  be  made  at  this  point  and  then  tracked  to  the  appropriate 

facets.  Because  it  is  the  calculation  of  the  scatter  type  and  process  that  is  the 

most  time  consuming  in  the  code,  this  would  lead  to  considerable  savings. 

* 

Another  variance  reduction  scheme  could  be  considered  when  the 
integration  is  being  performed  in  a  resonance  region.  Right  now,  the  draw  is 
uniform  over  the  incident  energy,  although  we  know  that  the  integral 
depends  on  In  a  group  with  a  large  number  of  resonances,  the  MC 

method  might  have  difficulty  in  converging  to  the  correct  answer.  Drawing 
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the  incident  energy  based  on  Gs(E')  might  improve  the  results  of  the 
integration. 

As  it  stands  the  code  runs  fairly  fast.  On  an  IBM  300MHz  Pentium  II 
processor  with  64Mbyte  RAM,  timing  runs  were  performed  for  a  single 
incident  group  and  facet.  For  the  21  x  40  facet  structure  (comparable  to  an 
enormous  Sae  level  symmetric  quadrature)  the  first  group  scatter  to  all 
scattered  groups  and  facets  was  calculated.  The  code  wall  clock  time  as  a 
function  of  draws  for  this  physics  model  is  plotted  in  Figure  50.  After  about 
one  hour,  80  million  T-matrix  samples  have  been  calculated. 


Run  Time  by  Number  of  Draws; 

(1  group,  1  facet)  ->  (30  groups,  762  facets) 


Figure  50.  Run  time  for  a  single  incident  group  and  facet. 

The  Sio  quadrature  consists  of  120  ordinates.  The  Cartesian  counterpart 
is  a  10  X  15  -^122  facet  quadrature.  As  will  be  shown  in  the  following 
chapter,  80  million  draws  are  not  required  for  accurate  transport  results. 


125 


Instead,  a  t5^ical  number  of  draws  for  a  full  set  of  data  is  about  five  million 
lasting  about  five  minutes.  That’s  5  incident  facets  by  30  groups  by  5 
minutes  each  or  a  complete  runtime  on  the  set  of  about  12.5  hours.  The  other 
advantage  is  that  this  calculation  is  only  performed  once  for  the  quadrature 
and  group  selected.  Stored  in  its  compacted  form  (say  on  a  CD),  rebalancing 
and  mixing  take  only  a  few  minutes  to  prepare  the  files  for  a  desired 
transport  problem. 

One  way  to  significantly  reduce  the  run  time  of  building  such  a  library  is 
to  run  each  incident  facet  and  group  calculation  on  a  separate  processor. 
Because  the  calculation  of  one  integration  is  completely  independent  of  the 
others,  the  data  can  be  combined  after  the  processors  have  completed  their 
calculations.  With  150  processors,  the  entire  122  facet  quadrature  library 
could  be  obtained  in  under  10  minutes  for  each  isotope. 

Finally  it  should  be  noted  that  the  real  test  of  how  well  these  libraries 
work  is  in  how  they  perform  in  the  transport  code.  Such  information  is 
presented  in  the  next  chapter.  The  following  section,  however,  continues  to 
look  at  the  scattering  T-matrix.  It  examines  how  the  MC  facet  method 
compares  with  its  spherical  harmonic  counterpart.  Particularly,  we 
demonstrate  the  non-negativity  of  the  MC  facet  method. 
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The  MC  Facet  Approach  vs.  the  Spherical  Harmonic  Method 

Up  to  this  point  we’ve  established  that  the  results  of  the  T-matrix  created 
with  the  MC  facet  method  are  consistent  with  our  expectation  of  the  physics. 
It  is  also  been  shown  that  the  MC  facet  method  can  be  implemented  to  create 
a  reasonably  accurate  T-matrix  in  a  modest  amount  of  time.  Such  a 
scattering  matrix  can  also  be  created  directly  from  the  spherical  harmonic 
expansion  of  the  scatter.  In  fact,  Tetran  performs  this  calculation  during  the 
transport  process.  We  can,  to  some  extent,  use  the  same  visualization  tools  to 
examine  the  behavior  of  the  spherical  harmonic  expansion  of  the  scattering 
cross  section. 

First,  examine  the  within  group  scatter  of  hydrogen  as  shown  in  Figure 
31.  The  P3,  P5  and  P7  Legendre  expansions  of  this  distribution  are  shown  in 
Figure  51. 
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Spherical  Harmonics  v.s.  T-Scat 


Figure  51  Legendre  Expansion  of  the  hydrogen  differential  scattering  cross 

section. 


The  effect  such  an  expansion  has  on  the  scattering  distribution  is  even 
more  dramatic  when  one  examines  the  spherical  harmonic.  Figure  52  shows 
the  P3  and  P7  spherical  harmonic  expansions  of  the  within  group  scatter 
(1-»1)  of  hydrogen  using  the  LANL-30  group  structure  (Preusser 
interpolation  of  the  scattered  data  with  TableCurve-3D).  The  chart  is  cut  off 
at  zero  to  show  exactly  where  the  truncated  expansion  goes  negative. 


128 


Hydrogen  P3  (1->1) 


Hydrogen  P7  (1->1) 


Figure  52.  Cartesian  map  of  the  P3  and  P7  spherical  harmonic  expansions  of 
iH.  Within  group  scatter  (1->^1)  LANL-30.  The  data  is  cut  off  at  zero  to 
demonstrate  the  sea  of  negativity. 

It  is  in  these  regions  that  the  potential  for  creating  negative  sources  exists. 
Increasing  the  order  of  the  polynomial  may  better  approximate  the  forward 
scattered  peak,  but  the  depth  of  its  negative  regions  becomes  more  severe  as 
was  shown  in  Figure  51. 

As  the  atomic  weight  of  the  target  nuclide  becomes  larger,  the  scattering 
probability  distribution  function,  may  have  more  structure. 

When  the  scatter  is  discrete  level  inelastic,  there  are  other  problems.  The 
third  order  polynomial  is  unable  to  match  the  structure  and  so  diffuses  the 
neutrons  in  all  directions.  Further,  as  expected,  this  approximation  has 
negative  components.  The  seventh  order  approximation  better  matches  the 
shape,  but  also  fails  to  maintain  non-negativity  (Figure  53). 
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i60(l->2)  P3 


i60(l->2)  P7 
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Figure  54.  Group  1^2  transfer  for  using  T-Scat  (some  interpolation  error 

in  74  facet  case). 

The  74  and  762  facet  cases  are  shown  for  comparison  in  Figure  54.  The 
762  facet  case  demonstrates  what  the  scattering  should  look  like.  The  74 
facet  case  (comparable  to  an  Ss  quadrature)  is  strictly  non-negative.  The 
source  facet  is  off  the  equator  at  ^  =  0.16  accounting  for  the  polar  asymmetry. 
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The  peaks  at  ^  «  0.75  are  a  bit  larger  than  the  actual  data,  and  are  artifacts 
of  the  interpolation  routine. 

Figure  55  shows  where  the  spherical  harmonic  approximations  go 
negative,  both  for  the  P3  and  P7  truncations.  The  bi-modal  behavior  of  the 
scatter  results  from  two  competing  mechanisms.  The  scatter  distribution, 
f{E[ is  highly  forward  peaked but  the  optimal  scattering  angle 

for  the  group-to-group  transfer  is  about  (see  Figure  38). 


Figure  55.  Spherical  harmonic  comparison  with  T-Scat  for  elastic  and 
inelastic  scatter  of  i^O.  (Group  l->2,  LiANL-30  Group  structure) 

Graphically,  it  is  easy  to  see  the  difficulties  with  using  Legendre 
polynomials  to  match  data  that  sometimes  is  only  C®  continuous.  It  is  a  more 
subtle  study  to  see  if  such  polynomials  adversely  effect  transport  codes.  The 
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abundance  of  attempts  to  fix  the  problem  with  source  negativity  suggests 
that  such  adverse  effects  exist. 

Researchers  at  AFIT  have  developed  a  new  spatial  characteristic  method 
that  allows  for  the  study  of  deep  penetration  problems  with  very  optically 
thick  cells^^-^®.  The  exponential  characteristic  method  (EC)  requires  that  the 
source  be  non-negative  or  the  code  will  faiF.  Other  transport  codes 
incorporating  exponential  methods  have  also  failed  when  dealing  with 
anisotropic  scattering^. 

Source  negativity  also  impacts  convergence  of  the  transport  codes.  A 
common  fix  is  to  set  the  source  to  zero  wherever  it  becomes  negative.  This 
also  impacts  the  convergence  of  the  code,  but  in  a  non-linear,  and  often 
unpredictable  way.  The  result  may  be  to  arrive  at  an  answer  that  is 
incorrect,  or  to  severely  impact  how  long  convergence  takes^.  The  effect  is 
more  dramatic  if  finer  group  structures  are  used.  In  the  examples  above,  a 
fairly  coarse  survey  structure  was  used  (LANL-30  Group). 

In  the  following  chapter  we  directly  compare  the  performance  of  the 
T-Scat  generated  T-matrix  with  its  spherical  harmonic  counterpart.  *-In  the 
next  section,  however,  we  continue  looking  at  the  scattering  matrix  generated 
by  T-Scat.  In  particular,  we  examine  the  problems  of  angular  support  for 
ordinate -to-ordinate  methods,  and  the  problems  of  numerical  diffusion  for  the 
facet  method. 
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The  Impact  of  the  Angular  Quadrature 

'  Earlier  we  mentioned  the  failure  of  ordinate-to-ordinate  methods  to 
scatter  to  the  next  group  if  the  nearest  scattering  angle  wasn’t  close  enough. 
We  quantify  that  here,  and  then  look  at  comparable  facet  quadratures.  We’ll 
find  that  the  MC  facet  method  guarantees  within  group  scatter  and  next 
group  scatter.  Further,  no  group-to- group  scatter  can  be  missed  if  it  is 
physically  allowed  to  happen. 

Recall  from  chapter  4  (Angular  Support)  that  the  scattering  angle  in  the 
laboratory  frame  is  given  by 


(188) 


For  elastic  scattering,  Q  =  0 ,  equation  (188)  can  be  re-written  as 


//i  E'-AE  ..  ■  E 

+  - (^-1) 


E' 


E'-AE 


(189) 


where  AE  is  the  energy  lost  in  the  colhsion.  Expanding  equation  (189)  in  a 
power  series  about  AE  =  0 ,  the  scattering  angle  for  small  energy  loss  is 


(190) 


Because  of  this  relation,  ordinate-to-ordinate  methods  lack  angular  support 
when  fine  energy  groups  are  used  with  modest  angular  quadratures. 
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Scattering  Within  Group 


We  first  examine  the  within  group  (1->1)  transfer.  Is  it  possible  to  scatter 
(change  directions)  while  remaining  in  the  same  group?  From  a  physical 
standpoint  the  answer  is  yes.  Neutrons  at  17  MeV  may  scatter  to  15.01  MeV 
and  still  remain  in  the  same  group  (group  1  of  the  LANL-SO  group  structure). 
Even  for  very  fine  energy  group  structures  a  loss  in  energy  involves  a  change 
in  direction.  But,  as  was  shown  earlier,  in  ordinate-to-ordinate  methods,  if 
there  is  no  directional  ordinate  within  the  scattering  range  (a  lack  of  angular 
support),  such  a  transfer  can  not  take  place. 

Using  equation  (31)  for  the  within  group  scatter  of  hydrogen, 

Ml  (1^1)  =  0.939.  To  guarantee  within  group  scatter  for  all  groups,  the  most 
stringent  requirement  is  for  group  2,  namely  ml  (2^2)  =  0.949.  The  next 

question  is  what  kind  of  angular  support  is  required  to  guarantee  within 
group  scatter? 

For  comparison,  we  use  a  level-symmetric  quadrature  to  represent  the 
ordinates  of  the  ordinate-to-ordinate  method.  Mikols  demonstrated  that 
Lobato  quadratures  are  even  more  stringent^h  The  nearest  neighbor  for  such 
a  quadrature  is  shown  in  Table  3  where  the  order  of  the  quadrature  and 
number  of  ordinates  in  the  quadrature  set  is  shown  next  to  the  nearest 
neighbor. 
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Table  3.  Nearest  Neighbor  for  Level  Symmetric  Quadratures. 


Level- 

Symmetric 

Sn 

#  of  Ordinates 

Nearest  jui^ 

S2 

4 

0.333 

S4 

24 

0.755 

Se 

48 

0.884 

Ss 

80 

0.922 

Sio 

120 

0.941 

Sl2 

168 

0.952 

Sl6 

288 

0.965 

To  achieve  within  group  scatter  for  all  groups  of  hydrogen  (using  the 
LANL-30  group  survey  structure)  the  order  must  be  at  least  12.  If  we 
perform  the  same  analysis  using  a  fine  group  structure  often  used  in 
transport  calculations  (e.g.,  Vitamin-J),  the  minimum  scattering  angle  is 
greater  than  0.9654  for  104  of  the  175  groups!  Because  increasing  the 
quadrature  order  indefinitely  becomes  untenable,  the  choice  must  be  made  to 
sacrifice  within  group  scatter,  ignoring  the  physics  of  the  calculation.  The 
MC  facet  method  has  no  such  restriction.  Because  the  minimum  angle 
between  adjacent  facets  is  0  within  group  scatter  is  guaranteed  for 

any  quadrature  set. 
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Angular  Support  for  Downscatter 

The  next  least  restrictive  case  is  to  examine  scatter  to  the  next  lower 
group.  How  close  must  the  angles  be  to  guarantee  scatter  to  the  next  lower 
group?  Again  we  examine  hydrogen.  The  maximum  energy  that  can  be  lost 
in  the  group -to- group  transfer  is  from  the  top  of  the  incident  group  to  the 
bottom  of  the  scattered  group.  For  the  LANL-30  group  structure,  (l->2) 
corresponds  to  17  MeV  ->  13.5  MeV.  The  result  is  a  scattering  cosine  of 

=  0.8912. 

Using  Table  3  again  we  see  that  an  Ss  quadrature  (80  ordinates)  is 
required.  For  the  Vitamin-J  structure,  15  of  the  175  angles  require  Sie  or 
higher.  Table  4  shows  the  number  of  groups  failing  as  the  quadrature  order 
is  reduced. 

Table  4.  Number  of  groups  Failing  Nearest  Neighbor  Downscatter  for 
Vitamin-J  group  structure  (175  total  groups). 


Level- 

Symmetric 

Sn 

#  of  Ordinates 

#  Failing  of  175 

S4 

24 

172 

Se 

48 

124 

Ss 

80 

109 

Sio 

120 

96 

Sl2 

168 

28 

Sl6 

288 

15 
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Of  course,  if  the  quadrature  fails  to  allow  downscattering  to  the  next 
group  when  the  physics  dictates  it  should,  that  group  is  simply  skipped.  The 
result  is  a  depressed  or  perhaps  non-existent  flux  in  the  lower  group. 

Again  the  MC  facet  method  guarantees  that  all  groups  are  scattered  to 
and  that  no  group  is  missed  because  the  scattering  is  continuous  in  angle 
within  the  calculation  of  the  T-matrix. 

As  an  illustration  of  this  t5q)e  of  facet  scatter,  consider  an  equal  area 
quadrature  as  previously  discussed.  Here,  the  scattering  is  strictly  from  the 
pole  to  the  facets  below.  We  first  examine  group  1^1  scatter  for  the  LANL- 
30  group  structure.  The  coarsest  structure  for  T-Scat  is  to  scatter  from  one 
hemisphere  to  the  next.  We  next  increase  the  number  of  total  latitudinal 
facets  to  2,  4,  10  and  100.  The  results  are  shown  in  Figure  56.  While  for  the 
fine  mesh  a  number  of  facets  are  reached,  even  the  coarsest  meshes,  as  is 
shown  by  equation  (35),  have  some  portion  of  the  scatter  going  to  the 
adjacent  facet. 
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Hydrogen  T-matrix  (1^1) 
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Hydrogen  T-matrix  (1->1) 
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Figure  56.  Guaranteed  within  group  scatter  [iH(l^l)]  for  the  MC  Facet 
Method.  Note  scale  change  in  100  latitude  chart. 

As  expected  the  MC  facet  guarantees  scatter  to  groups  where  the  physics 
says  they  should.  Again,  this  is  because  at  the  heart  of  the  MC  facet  method 
the  scattering  is  taken  as  continuous  in  both  energy  and  angle. 


Skipping  Groups 

Finally,  it  should  be  noted  that  ordinate-to-ordinate  methods  may  also 
skip  groups  if  there  is  insufficient  angular  support.  Because  the  MC  facet 
method  guarantees  scatter  to  all  e  [-I,  l] ,  no  groups  will  be  skipped  if  the 

physics  dictates  they’re  to  be  scattered  to. 
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Angular  Dispersion 


The  price  paid  for  guaranteeing  within  group,  next  group  and  non-skipped 
group  scatter  is  in  angular  dispersion.  Figure  57  shows  such  dispersion  by 
comparing  facet  structures  (the  effect  is  also  noticeable  in  Figure  56).  Both 
are  for  hydrogen  within  group  scatter.  The  coarsest  structure  tends  to 
redistribute  a  small  portion  of  the  scatter  into  a  direction  that  by  kinematics 
cannot  be  reached  (recall  that  for  the  polecap  quadrature  there  are  ordinates 
at  1  and  -1,  while  the  other  ordinates  are  centered  on  the  facet). 


Facet  Angular  Dispersion 


Figure  57.  Facet  angular  dispersion  for  ^H(l^l)  scatter  (expanded  axis). 

Of  course,  Legendre  polynomials  suffer  from  the  same  type  of  dispersion  by 
assuming  scatter  is  possible  in  regions  where  it  is  not.  This  is  evident  by 
comparing  the  P3  and  P7  approximations  in  Figure  51.  Further,  such 
approximations  must  have  negative  regions  where  the  actual  matrix  is  zero. 
Just  as  the  facet  method  is  dependent  on  the  quadrature  used,  the  spherical 
harmonic  expansion  is  also  dependent  on  the  order  of  the  Sn  method  used. 
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The  best  way  to  determine  if  such  dispersion  is  problematic  is  to  try  the 
cross-section  data  libraries  in  a  transport  code.  This  is  done  in  the  following 
chapter  which  evaluates  the  transport  code  results.  In  particular  we’ll 
examine  comparable  quadratures  to  determine  if  there’s  a  significant 
problem  or  not. 

One  way  to  reduce  the  angular  dispersion  in  the  facet  method  is  to 
perform  an  angular  rebalance  to  conserve  the  first  angular  moment  of  the 
scatter.  This  is  discussed  in  the  next  section. 

The  Effects  of  Angular  Rebalance 

Ideally  we  would  examine  the  effects  of  angular  rebalance  on  diffusion- 
like  problems  in  the  transport  code  chapter.  Unfortunately,  the  absence  of  a 
convergence  accelerator  in  TETRAN  makes  such  an  examination  impossible 
at  the  present  time.  In  this  section  we’ll  continue  to  examine  the  T-matrix 
using  the  tools  developed  earlier.  In  particular,  we  want  to  see  that  the 
rebalance  doesn’t  turn  the  matrix  into  something  unrealistic.  And,  better,  we 
hope  that  the  effects  of  the  rebalance  are  consistent  with  our  expectations. 

Perhaps  the  best  benefit  of  rebalance  is  the  reduction  of  angular 
dispersion  as  illustrated  in  Figure  57.  In  that  figure,  the  T-matrix  values 
before  and  after  rebalancing  are  shown.  In  general,  we  would  hope  that, 
after  the  rebalance,  the  peaks  of  the  scattering  distributions  would  be  more 
pronounced  and  the  dispersive  effects  reduced. 
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To  look  at  this,  take  a  standard  74  facet  quadrature  set  with  incident 
neutrons  through  the  polecap.  The  azimuthal  angle  is  collapsed  and  the 
T-matrix  results  for  before  and  after  rebalancing  are  shown  in  Figure  58. 
The  purpose  of  rebalancing,  as  described  earlier,  is  to  conserve  the  first 
angular  moment  of  the  scatter.  This  usually  results  in  redistributing  the 
scatter  matrices  towards  peaks. 


Figure  58.  Reduction  of  dispersion  using  angular  rebalancing.  Here,  74  facet 

quadrature  for  (1-»1) 

The  effect  can  be  seen  more  dramatically  in  three  dimensions  by 
examining  the  Cartesian  grid  of  similar  data.  Here  a  9  x  12  ->  86  facet 
structure  is  used  to  center  the  source  of  the  scatter  at  ^  =  0  (Figure  59). 
Again,  the  T-matrix  values  before  and  after  rebalancing  are  presented  side  by 
side  for  hydrogen.  For  the  within  group  scatter  (1-^1),  we  see  that  the  peak 
of  the  rebalanced  data  is  higher  and  better  defined. 
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Figure  59.  Effects  of  angular  rebalance  for  hydrogen  (1->1). 


Figure  60.  Effect  of  angular  rebalance  for  equatorial  facets,  (1->1). 


Next  we  examine  the  same  type  of  results  for  the  oxygen  (l->2)  scatter. 
As  for  hydrogen,  we  look  at  a  linear  plot  of  the  zeroth  moments  of  the  facets 
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along  the  equator  before  and  after  rebalancing.  Rebalance  effects  are 
compared  to  a  baseline  of  100  equatorial  facets  that  show  the  strong 
structure  of  the  ^^0  (l->2)  scatter.  Rebalancing  of  the  coarse  mesh  (9xl2->86 
facet  quadrature)  and  the  medium  mesh  (9x20^142  facet  quadrature)  are 
shown  in  Figure  61  and  Figure  62  respectively. 

Conservation  of  the  first  angular  moment  reduced  the  structure  in  the 
rebalanced  coarse  mesh.  For  the  medium  mesh,  just  as  it  was  for  the  within 
group  scatter  of  hydrogen,  the  peaks  are  enhanced,  more  closely  resembling 
the  actual  scatter  structure. 


Figure  61.  Rebalance  of  (l->2)  T-matrix  for  12  equatorial  facets. 
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Figure  62.  Rebalance  of  (l->2)  T-matrix  for  20  equatorial  facets. 


Overall,  the  rebalance  seems  to  reduce  dispersion,  particularly  for  modest 
mesh  sizes.  This  is,  however,  the  secondary  goal  of  the  rebalance.  Its 
primary  mission  is  to  ensure  that  the  scattering  libraries  perform  well  for 
diffusion-like  problems.  Our  current  code  limitations  (TETRAN’s  lack  of  a 
convergence  accelerator)  do  not  allow  us  to  examine  the  libraries  at  low 
energies,  but  the  results  shown  here  suggest  that  rebalancing  will,  at  the 
very  least,  not  impair  the  accuracy  and  may  even  improve  it.  Even  for  high 
energy  problems  well  see  that  rebalancing  can,  in  fact,  improve  the  accuracy 
of  our  transport  code  results. 


Scattering  Matrix  Conclusions 

After  a  thorough  examination  of  the  scattering  matrix  weVe  shown  a 
number  of  things:  The  MC  facet  method  creates  a  T-matrix  that  is 
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•  consistent  with  the  physics, 

•  guaranteed  non-negative, 

•  combines  both  elastic  and  inelastic  collisions, 

•  guarantees  within  group  and  next  group  scatter  regardless  of 
quadrature  or  energy  group  structure,  and 

•  can  be  rebalanced  to  conserve  the  first  angular  moment  of  the  scatter. 

Further,  it  does  so  at  a  reasonable  accuracy  and  in  a  reasonable  amount  of 

time.  And,  if  the  user  is  willing  to  spend  the  time,  noting  that  the  library 
need  only  be  built  once,  any  T-matrix  accuracy  can  be  achieved. 

While  significant  achievements  on  their  own,  it  remains  to  be  seen  how 
well  such  a  T-matrix  will  work  in  a  transport  code.  This  is  the  subject  of  the 
following  chapter.  And  although  the  Cartesian  angular  quadrature  is  not 
ideal,  the  next  chapter  will  provide  some  comparison  between  the  facet 
method  and  the  spherical  harmonic  approach  and  demonstrate  the  viabihty 
of  the  MC  facet  method  for  neutral  particle  transport. 
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10.  Transport  Code  Performance 


So  far,  we’ve  discussed  the  theory  for  the  MC  facet  method  and  the  code 
that  implements  that  theory  to  create  a  T-matrix.  Examination  of  this 
T-matrix  data  shows  that  it  models  the  physics,  is  non-negative,  guarantees 
within-group  and  next-group  scatter,  and  conserves  the  first  angular  moment 
of  the  scatter.  This  chapter  now  looks  at  how  well  such  a  cross-section  library 
performs  in  a  transport  code.  We’ll  find  that  the  T-matrix  created  with  the 
MC  facet  method  performs  well  for  a  number  of  sample  problems,  and  in 
most  cases  is  superior  to  a  comparable  spherical  harmonic  cross-section 
library. 

This  chapter  begins  to  explore  the  capabilities  of  the  T-matrix  approach. 
Our  goal  is  to  demonstrate  the  capabilities  of  the  T-matrices  created  by 
T-Scat  in  the  transport  code.  These  include 

•  Non-negative  T-matrix  libraries  for  single  elements  (e.g.,  ^H,  i®0). 

•  Non-negative  T-matrix  libraries  for  compounds  (e.g.,  H2O,  BO2) 

•  Support  for  geometries  with  multiple  materials  (e.g.,  H2O,  Fe,...). 

Also  investigated  is  the  impact  of  various  quadrature  and  group  structures 
for  some  of  the  problems  described  and  comparisons  of  comparable  changes  to 
the  spherical  harmonic  counterpart. 
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Overview  of  the  Transport  Codes 


The  primary  code  used  to  investigate  the  new  non-negative  cross  section 
libraries  was  TETRAN.  TETRAN  was  developed  at  AFIT^.  It  is  an 
unstructured  tetrahedra  mesh  discrete  ordinates  radiation  transport  code. 
TETRAN  solves  the  linear  time-independent  Boltzmann  transport  equation 
using  linear  or  exponential  characteristic  spatial  quadratures  and  is  capable 
of  working  with  multi-group  cross  sections  and  general  anisotropic  scattering. 

As  discussed  earlier,  TETRAN  creates  a  T-matrix  using  spherical 
harmonics.  We  modified  the  code  so  that  it  could  read  T-matrix  data  directly 
from  the  cross  section  libraries  created  by  T-Scat,  T-Balancer  and  T-Mixer. 

In  this  way  we  can  directly  compare  the  results  of  the  two  methods  using  the 
same  spatial  mesh  and  group  structure.  Of  course  the  angular  quadratures 
are  different,  but  effort  is  made  to  keep  such  quadratures  as  similar  as 
possible  (e.g.,  same  number  of  ordinates). 

Data  is  also  compared  with  the  results  of  MCNP  using  continuous  energy 
and  angle  ACE  libraries.  These  are  the  same  libraries  used  to  calculate  the 
cross  section  data  within  T-Scat. 

On  occasion  we  also  introduce  comparison  results  from  MCSN.  As 
discussed  earlier,  this  code  performs  transport  continuous  in  space  but 
discrete  in  energy  and  angle.  It  too  reads  the  MC  facet  method  T-matrix 
directly  from  file  to  create  its  cumulative  distribution  functions,  for  the  MC 
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transport  process.  MCSN,  TETRAN,  and  MCNP  provide  an  excellent  set  of 
tools  to  evaluate  how  well  this  new  method  performs. 

Nuclides  Investigated 

T-Scat  currently  implements  elastic,  level  inelastic  and  evaporation 
spectrum  scattering  laws.  Nuclides  and  energy  ranges  were  selected  so  that 
such  scattering  laws  were  the  only  ones  required.  Nuclides  were  also  selected 
to  demonstrate  the  anisotropic  capabilities  of  the  new  method.  Hence,  lighter 
nuclides  were  of  greatest  interest.  Nuclides  and  mixes  tested  include  ^H,  “^Li, 
lOR,  16  O,  27A1, 56Fe,  H2O,  BO2,  AI2O3,  and  others.  The  results  presented  in 
this  chapter  are  representative  of  the  results  found  for  all  these  materials. 
Where  possible,  to  demonstrate  operation  at  high  anisotropy,  the  highest 
energy  groups  of  the  LANL-30  and  VITAMIN-J  structures  are  used.  In  some 
instances  lower  energy  groups  are  used  to  avoid  unavailable  types  of  scatter 
[e.g.,  (n,  2n)]. 

Definition  of  Terms 
Units 

The  output  data  of  TETRAN  is  consistent  with  the  tally  output  of  MCNP. 
Such  tallies  are  normalized  to  be  per  source  particle.  TETRAN  provides 
results  for  the 

1.  current  through  a  surface  (particles). 


148 


2.  flux  averaged  over  a  surface  (particles/cm^),  and 

3.  flux  averaged  over  a  material  region  (particles/cm^). 

If  the  source  is  given  in  particles  per  unit  time,  so  too  are  the  current  and 
fluxes  described  above.  In  this  research,  only  box  geometries  are  used.  All 
external  surfaces  are  surrounded  by  vacuum.  A  cube  of  water  represents  a 
single  material  region.  The  water-iron-water  geometry  used  in  this  research 
is  a  cube  in  a  cube  in  a  cube  and  represents  three  regions. 

Error 

MCNP  reports  the  relative  error  for  each  of  its  estimates  of  the  mean. 

The  absolute  error  is  based  on  the  variance  of  the  estimated  mean  calculated 
during  the  tally.  The  discrete  ordinates  code  TETRAN  converges  to  a 
tolerance  level  of  10  ®.  When  an  error  is  shown  for  MCNP,  it  is  the  relative 
error  reported  by  the  code,  or  absolute  error  defined  as  the  relative  error 
times  the  mean.  Errors  reported  for  TETRAN  and  MCSN  are  based  on  the 
difference  between  their  calculated  values  and  those  of  MCNP.  As  mentioned 
earlier,  it  is  unreasonable  to  expect  that  such  codes  would  converge  to  the 
exact  same  result.  What  we  hope  to  see  is  that  MCSN  lies  between  the 
results  of  TETRAN  and  MCNP.  It  is  also  hoped  that  results  using  the  MC 
facet  approach  perform  somewhat  closer  to  the  MCNP  result  than  those  using 
spherical  harmonics.  It  will  be  shown  that,  in  most  instances,  this  is  indeed 
the  case. 
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Spatial  Quadratures 

Of  course  neither  MCNP  nor  MCSN  use  spatial  quadratures.  TETRAN 
may  implement  one  of  two:  linear  characteristic  (LC)  and  exponential 
characteristic  (EC). 

The  linear  characteristic  method,  as  developed  by  Mathews®®,  assumes  a 
linear  source  and  inflow  flux  distribution.  Spatial  moments  are  taken  of  the 
LC  source  distribution  and  a  system  of  equations  is  set  up  to  solve  for  the 
source.  The  LC  method  performs  well  if  the  spatial  cells  are  optically  not  too 
thick.  The  method,  however,  may  return  negative  flux  values  even  if  the 
cross  section  data  is  strictly  positive. 

The  exponential  characteristic  method,  however,  creates  strictly  non¬ 
negative  fluxes.  The  EC  method  takes  zeroth  and  first  spatial  moments  of 
the  characteristic  solution  to  the  BTE  and  an  assumed  exponential 
distribution  of  the  scattering  source  in  a  cell  to  calculate  average  and  first 
moments  of  the  angular  flux®®.  The  method  is  excellent  for  deep  penetration 
problems.  The  method,  however,  requires  the  source  be  non-negative.  The 
availability  of  an  anisotropic,  non-negative,  cross  section  library  (the*- thrust  of 
this  research)  makes  the  implementation  of  the  EC  method  now  feasible  for 
multi-group,  anisotropic  problems.  Both  these  methods  are  discussed  in 
detail  in  Miller’s  dissertation'^. 

Because  transport  calculations  implementing  the  LC  method  are  quicker, 
many  of  the  results  of  this  chapter  use  LC.  Use  of  the  EC  method  is 
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demonstrated  for  some  problems  and  in  particular  for  problems  with  optically 
thick  cells  where  the  LC  method  fails  to  provide  positive  fluxes. 


Uniform  Cube  Test  Problems 


The  first  test  of  the  T-Scat  data  was  with  a  uniform  cube.  The  cube 
consists  of  a  single  material  (10cm  x  10cm  x  10cm).  A  source  emits 
isotropically  throughout  the  material  (1  particle/cm^-sec)  and  uniformly 
within  a  single  energy  group  of  the  given  group  structure. 


Figure  63.  Dimensions  of  the  simple  cube. 


Introducing  Data  for  Single  Nuclides 

This  section  limits  the  cube  material  to  a  single  nuclide.  It  introduces 
some  of  the  results  from  TETRAN  and  MCNP  and  ends  with  a  presentation 
of  data  from  MCSN.  In  so  doing,  a  number  of  characteristics  of  the  MC  facet 
method  are  revealed.  Most  comparisons  are  based  on  the  region  average 
scalar  flux 


JlOcm  flOcm  rlOcm  ,  /  \ 

0  dz<l>g{x,y,z)  . 


lOOOcm^ 


(191) 
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The  first  material  tested  was  free  hydrogen  (1.0  g/cm^)  at  300K.  The 
problem  was  calculated  using  MCNP,  TETRAN  with  spherical 
harmonics  (SH),  and  TETRAN  with  the  T-Scat  library  (five  milhon  scatter 
draws).  This  density  is,  of  course,  unrealistic.  It  was  used  for  convenience, 
although  one  could  use  a  larger  cube  at  lower  density  to  obtain  the  same 
results.  The  cube  is  divided  into  161  spatial  cells  (tetrahedra).  The  LANL-30 
group  structure  was  used  and  the  group  scalar  flux  averaged  over  the 
material  (equation  (191))  for  the  first  10  groups  is  shown  in  Figure  64.  The 
total  group  scalar  flux  spectrum  is  found  by  dividing  the  group  scalar  flux  by 
the  width  of  the  group  in  MeV 

(192) 

This  is  shown  in  Figure  65.  The  SH  method  employs  an  S8  quadrature 
(80  ordinates)  with  a  P7  approximation  while  T-Scat  uses  the  8  x  12  ->  74 
facet  structure.  On  this  scale  it  is  difficult  to  see  exactly  which  of  the 
TETRAN  results  is  closer  to  the  MCNP  solution.  What  can  be  said  is  that  the 
T-Scat  hbrary  is  performing  well.  We  can  more  closely  examine  the  error  of 
these  data  by  comparing  them  to  the  MCNP  result. 
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Figure  64.  Hydrogen  cube  test  problem  comparing  volume  scalar  flux. 

Again,  as  noted  above,  the  uncertainty  reported  for  MCNP  is  its  own 
reported  relative  uncertainty.  We  also  provide  the  99%  confidence  level  of 
the  MCNP  uncertainty  (fggo/^  =  2.58s)  for  comparison.  The  differences 

reported  for  TETRAN  are  based  on  their  relative  differences  to  the  MCNP 
value  namely 

1 AA  W  TETRAN  ~  ^ MCNP  n/  /I  Q  0\ 

S  TETRAN  =  100^ - - - '-%  (193) 

T  MCNP 

The  purpose  of  reporting  these  values  is  to  give  the  reader  confidence  that 
the  MC  facet  method  performs  reasonably  well  in  a  transport  code.  The 
method  is  certainly  not  optimized  at  this  time  (particularly  in  its  angular 
quadrature).  The  differences  should  act  as  a  guide  in  assessing  the  feasibility 
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of  the  MC  facet  method  noting  that  tjT)ical  cross  section  library  uncertainties 
can  be  on  the  order  of  a  few  percent. 


17  15  13.5  12  10  7.79  6.07  3.68  2.865  2.232 

Group  Max  (MeV) 


Figure  65.  Hydrogen  cube  volume  scalar  flux  spectrum. 


Figure  66.  Relative  difference  between  TETRAN  methods  and  MCNP  (^H). 


Changes  from  Angular  Quadrature 

Figure  66  charts  the  relative  difference  of  the  various  methods.  Even  with 
six  fewer  ordinates,  the  T-matrix  library  calculated  with  T-Scat  has  a 
generally  lower  difference  than  its  spherical  harmonic  counterpart.  The  lines 
in  the  chart  are  added  only  for  clarity  and  do  not  represent  any  data. 

The  transport  code  performance  is  improved  if  the  number  of  directional 
ordinates  using  T-Scat  is  increased  from  six  fewer  than  the  S8  to  six  more. 
Using  a  8  x  14  ->  86  facet  structure  the  new  values  are  shown  in  Figure  67. 


Figure  67.  Reducing  the  T-Scat  difference  by  increasing  from  74  to  86  facets. 

We  continue  to  examine  angular  quadratures  sets  with  more  ordinates  by 
performing  the  same  type  of  calculation  for  The  SH  was  an  S10/P3  and 
T-Scat  used  114  facets.  The  results  are  shown  in  Figure  68  and  demonstrate 
again  that  T-Scat  is  providing  quality  results. 
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Figure  68.  Simple  cube  for  shows  marked  improvement  for  T-Scat. 


Figure  69.  Simple  cube  for  shows  marked  improvement  for  T-Scat. 


Group  Structure  Effects 

We  next  examined  the  hydrogen  cube  using  the  finer  group  structure 
VITAMIN- J.  While  both  methods  under  predict  the  scatter  to  the  next  group, 
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%  Difference 


Of  course  reducing  the  angular  quadrature  can  also  impact  the  error.  We 


next  compare  an  SH  S4/P7  with  a  6  x  6  26  facet  quadrature. 


Scalar  Flux  Spectrum  for  Hydrogen  Cube 
(Low  number  of  directional  ordinates) 


Group  Max  (MeV) 


Figure  72.  Scalar  flux  for  a  small  number  of  directional  ordinates. 


Figure  73.  Flux  differences  for  a  small  number  of  directional  ordinates. 

Figure  72  shows  the  average  scalar  flux  in  the  cube  while  Figure  73  shows 
the  differences.  Of  particular  interest  is  the  improved  performance  of  the 


158 


facet  method  in  the  group  below  the  source  group.  While  both  methods 
under  predict  the  value  of  the  flux  in  the  second  group,  the  facet  method  does 
provide  about  a  factor  of  two  improvement. 

Spatial  Mesh  Variations 

Next,  we  look  at  what  happens  when  the  spatial  mesh  is  reduced.  As  for 
the  angular  quadrature  reduction,  this  is  done  simply  to  verify  appropriate 
operation  of  the  scattering  library.  Increasing  the  number  of  tetrahedra 
reduces  the  relative  differences  with  respect  to  the  MCNP  calculation.  Such 
reductions,  however,  are  negligible  suggesting  that  the  coarse  mesh  (161 
cells)  used  previously  was  sufficiently  fine). 


Figure  74.  Error  reduction  using  finer  spatial  mesh. 
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Introducing  the  MCSN  Method 


We  next  introduce  the  results  from  MCSN.  This  is  a  new  code  developed 
at  AFIT25,  It  is  spatially  continuous  but  discrete  in  angle  and  energy.  As 
mentioned  earlier,  it  requires  a  positive  scattering  matrix  to  build  the 
cumulative  distributions  used  to  sample  the  scattered  energy  and  facet. 
Because  it  uses  the  exact  same  T-matrix  library  as  TETRAN  (that  created  by 
the  MC  facet  method),  it  is  expected  that  the  commonality  of  being 
continuous  in  space  would  place  its  answers  closer  to  those  of  MCNP  than 
TETRAN.  In  fact,  MCSN  should  give  the  converged  result  that  TETRAN 
would  approach  as  the  spatial  mesh  is  refined.  Like  MCNP,  MCSN  reports 
its  own  relative  uncertainty.  For  the  calculations  presented  here,  such 
uncertainty  is  about  an  order  of  magnitude  less  than  the  comparable 
uncertainty  of  MCNP.  Figure  75  shows  that  the  differences  between 
TETRAN  and  MCSN  are  slight.  A  look  at  the  relative  differences  reveals 
that  our  predictions  were  correct.  The  MCSN  calculation  lies  somewhat 
closer  to  the  MCNP  result  than  the  TETRAN  calculation. 
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Figure  75.  MCSN  applied  to  the  161  cell  hydrogen  cube. 


Figure  76.  Error  of  TETRAN  and  MCSN  methods  using  T-Scat  library. 

We  can  also  examine  the  spatial  mesh  in  a  plot  similar  to  Figure  74. 
Here,  the  relative  difference  between  MCSN  and  TETRAN  is  given  where 
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both  methods  use  the  same  angular  quadrature  (74  facet)  but  TETRAN  is  run 
using  two  different  spatial  meshes  as  before  (161  and  1292  cells).  The  data  is 
plotted  along  with  the  MCNP  result  in  Figure  77  , 


Spatial  Convergence  using  MCSN 
74  Facet  Quadrature 


Group  Max  (MeV) 

Figure  77.  Difference  between  MCSN  and  TETRAN  for  two  spatial  meshes. 

It  is  not  the  purpose  of  this  research  to  demonstrate  the  capabilities  of  the 
MCSN  method.  Rather,  it  is  to  show  that  the  MC  facet  method  has  made  it 
possible  for  such  a  code  to  be  implemented,  and  that  its  results  are  consistent 
with  expectations.  Even  with  the  quadrature  limitations  imposed  on  the 
T-matrix  library,  results  are  excellent  for  the  average  scalar  flux  calculations. 
Similar  cube  calculations  were  performed  using  other  nuclides,  C^Li,  ^^Al,  ^^O) 
with  similar  results.  We  examine  a  number  of  mixed  materials  next,  while 
exploring  other  aspects  of  the  T-Scat  libraries. 


162 


Mixed  Nuclides  in  a  Single  Material 


Does  the  number  of  draws  taken  to  calculate  the  T-matrix  impact  the 
results  of  the  transport  code?  The  answer  is  yes.  The  same  cube  geometry 
above  was  used.  This  time  the  material  was  Li20  (1  g/cm^).  Again  the 
material  contained  an  isotropic  source  of  1  neutron/cm^-s.  The  T-matrix 
quadrature  was  74  facets,  and  the  LANL-30  group  structure  was  used.  The 
matrices  were  built  with  320  thousand,  1.28  million  and  5.12  million  draws 
and  then  run  in  TETRAN.  The  draws  were  independent,  each  started  with 
different  seeds  for  the  random  number  generator.  The  error  data  is  compared 
in  Figure  78. 


Flux  Differences  byT-Scat  Draws 
Li20  Cube 
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Figure  78.  Transport  error  may  depend  on  the  number  of  draws  to  calculate 

the  T-matrix  library. 


163 


For  the  simple  cube,  the  flux  differences  are  inconsequential.  A  closer 
look  at  the  second  group’s  relation  to  the  MCNP  solution  is  shown  in 
Figure  79. 


2nd  Group  Comparison  with  MCNP 


Figure  79.  Close  look  at  group  two  relation  to  MCNP  as  the  number  of  draws 
is  increased  in  creating  the  scattering  library. 


The  relative  error  of  the  MCNP  solution  was  0.05%  and  is  shown  by  the 
error  bars  in  the  figure.  Typically,  we  have  observed  little  practical 
improvement  beyond  5  million  draws  for  each  incident  facet  in  each  incident 
group.  This  is  true  for  the  scattering  matrices  used  in  all  of  the  geometries 
and  mixes  tested.  Unless  otherwise  stated,  it  is  the  default  number  of  draws 
taken  to  calculate  the  T-Matrix  as  described  earlier. 

The  next  mixture  of  nuclides  examined  is  the  most  common  --  water. 
Figure  80  shows  the  flux  calculation  for  water  using  the  four  transport 
methods  discussed  earlier.  The  error  calculation  reveals  results  consistent 
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with  what  was  shown  earlier.  The  MCSN  results  have  a  little  less  error  than 
TETRAN  results  using  T-Scat,  which  in-turn  is  less  than  the  errors  from  the 
SH  S8/P3  calculation. 


Figure  80.  Flux  comparisons  for  the  water  cube. 


Figure  81.  Error  of  the  methods  in  the  water  cube. 
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As  mentioned  earlier,  MCSN,  like  MCNP,  reports  its  own  relative 
uncertainty.  Figure  82  shows  how  the  two  compare  for  the  previous  example. 


Figure  82.  Stochastic  errors  of  MCSN  and  MCNP  in  the  H2O  cube. 


Rebalance  Effects  in  the  Transport  Process 

As  previously  discussed,  it  is  expected  that  rebalancing  the  first  angular 
moment  of  the  scatter  would  have  little  impact  on  high  energy  problems.  The 
rebalance  is  meant  primarily  for  diffusion  like  problems.  Figure  83  shows 
that  this  is  indeed  the  case.  Angular  rebalance  for  this  problem  made 
essentially  no  substantive  changes  to  the  scaler  flux  calculations  of  the 
transport  code. 
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Flux  Error  for  Hydrogen  Cube 


17  15  13.5  12  10  7.79  6.07  3.68  2.865  2.232 

Group  M  ax  (M  eV) 

Figure  83.  Effects  of  rebalancing  74  facets  for  the  hydrogen  cube. 

Currents  and  Quadratures 

The  level  symmetric  quadrature  used  here  for  the  spherical  harmonic 
expansions  provides  no  preferential  treatment  for  any  axis.  The  ordinate 
directions  are  invariant  to  90°  rotations  about  any  axis.  Current 
measurements  on  cube  faces  should,  therefore,  have  essentially  the  same 
values.  Unfortunately,  the  spatial  mesh  used  for  the  TETRAN  code  does  not 
have  this  same  property.  There  does  exist  some  variation  of  currents  on  the 
cube  faces  because  of  the  variation  in  the  spatial  mesh.  The  Cartesian 
quadrature  (with  pole  caps)  used  for  the  facet  quadratures  also  has  no  such 
symmetry.  It  is  invariant  to  rotations  about  the  s-axis  only,  and  to  reflections 
across  the  x-y  plane.  Hence,  for  a  cube  with  no  spatial  mesh  asymmetries. 


167 


currents  should  be  the  same  for  the  x'*'  ,x  ,y'^,  and  y  faces  of  a  cube,  but 

may  be  different  from  the  z'^  and  faces.  The  partial  current  in  the  +x 
direction  is 


Jx=  (194) 

with  the  unit  normal  in  the  positive  x  direction.  The  same  equation  holds 

for  the  other  directions.  For  an  isotropic  source  near  the  center  of  the  cube 
and  equal  quadrature  weights  (wn)  the  only  variabihty  is  in  the  magnitude  of 

^4  facet  quadrature, 

E  Z  (I,  n„)=i8.7.  (196) 

[e„A.>0]  [e^n,i>0] 


and  ^  •Q„)=18.8.  (196) 

This  suggests  a  bias  in  the  quadrature  that  will  increase  the  z  currents  such 
that,  for  the  cube  geometry 

j^^=j;>j^  =  j;=j;=j-  .  (197) 

With  such  a  bias,  we  expect  the  integral  current  out  the  ±z  faces  of  the 
cube  to  be  larger  than  the  integral  current  out  the  ±x  and  ±y  faces  of  the  cube. 
These  integral  values  are  reported  in  TETRAN’s  summary  file.  The  results 
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for  the  161  tetrahedron  water  cube  are  presented  in  Figure  84  and  Figure  85. 


Here,  the  1®*  and  2“^  group  currents  are  shown. 


Group  1  Surface  Currents:  H2O 


Figure  84.  Total  group  1  cube  surface  currents  of  H2O  cube  show  asymmetry 

in  facet  quadrature. 


Figure  85.  Total  ^oup  2  cube  surface  currents  of  H2O  cube  show  asymmetry 

in  facet  quadrature. 
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The  figures  show  that  level  symmetric  quadrature  does  produce  relatively 
uniform  currents  through  the  six  faces  suggesting  the  primary  variation  of 
the  T-Scat  method  is  due  to  the  angular  quadrature.  The  facet  quadrature, 
however,  results  in  a  predictable  difference  between  the  azimuthal  faces 
(ac  and  y)  and  polar  (z)  faces.  For  the  most  part  the  T-Scat  data  gives  better 
agreement  to  the  MCNP  calculation. 

To  investigate  other  aspects  of  the  MC  facet  method  for  anisotropic  scatter 
and  its  implications  for  use  in  sample  transport  problems,  a  more 
complicated  geometry  is  needed.  Beyond  demonstrating  the  efficacy  of  the 
method,  we’ll  also  demonstrate,  for  the  first  time,  the  use  of  the  EC  transport 
using  multi-group  anisotropic  scattering  matrices. 

Nested  Cube  Test  Problems 

These  problems  are  based  on  a  test  problem  used  by  Miller'^.  The 
geometry  is  three  cubes,  nested  centrally  as  shown  in  Figure  86.  There  are 
three  regions: 

A.  Source  region— within  innermost  cube  (material  1). 

B.  Middle  region-annulus  between  innermost  and  middle  cube 
(material  2). 

C.  Outer  region-annulus  between  middle  and  outer  cubes  (material  1). 

In  this  section  we  present  results  for  two  material  combinations:  Water- 

Iron-Water  and  Water-Boron-Water.  In  so  doing,  we  will  demonstrate  the 
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use  of  multiple  materials,  developed  by  T-Scat,  for  use  in  the  same  transport 
problem.  We  include  the  use  of  a  library  consisting  of  mixed  nuclides  for 

water  using  \H  and  ^|0  data.  Finally,  we  also  examine  the  effects  of  deep 
penetration  and  have  the  opportunity  to  demonstrate  (for  the  first  time)  the 
EC  method  using  anisotropic  scatter. 


Figure  86.  Mesh  and  dimensions  of  the  cube-cube-cube  geometry. 

We  first  examine  the  water-iron- water  problem.  For  simplicity,  H2O  (as 
described  above)  and  the  single  nuclide  are  the  materials  used. 
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Figure  87.  Water-iron-water  fluxes  for  six  groups. 


The  effects  of  the  penetration  through  the  shield  are  shown  for  each  group 
in  sequence.  To  avoid  unsupported  scatters,  the  source  emits  uniformly  in 
group  8  of  the  LANL-SO  group  structure  (2.865-3.68  MeV).  The  various 
methods  (MCNP,  TETRAN,  etc. )  are  placed  side  by  side  for  direct 
comparison.  Differences  with  MCNP,  as  before,  are  shown  as  line  charts. 
Figure  88  gives  a  quick  look  at  the  differences  for  this  problem  using  the 
T-Scat  library,  and  TETRAN’s  LC  spatial  quadrature  on  528  spatial  cells  and 
74  facets  (or  Sg).  Detailed  comparisons  with  MCNP  for  each  region  are 
shown  later. 
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Differences  of  the  Scalar  Flux  by  Group 
T-Scat  (LC) 


Figure  88  Presentation  of  H20-Fe-H20  differences  with  MCNP  by  region. 

The  optical  thickness  of  the  cells  plays  an  important  role  for  this  problem. 
The  optical  thickness  of  a  cell  is 

s  =  cr^cm”^  j  Z(cm)  (198) 

where  cr  is  the  macroscopic  total  cross-section,  and  I  is  the  distance  a  neutron 
travels  to  cross  the  tetrahedron.  The  exponential  characteristic  method  was 
designed  to  handle  optically  thick  cells.  The  average  optical  thicknesses  for 
the  tetrahedra  in  this  problem  are  shown  in  Table  5. 


Table  5.  Average  cell  optical  thickness  for  H20-®®Fe-H20  system. 


Group  (MeV) 

Water/Source 

®®Fe  (7.86g/cm3) 

Water/Outer 

3.68  -  2.865 

1.96 

3.233 

5.016 

2.865  -  2.232 

1.75 

3.173 

4.482 

2.232  -  1.738 

2.238 

2.935 

5.723 

1.738  -  1.353 

2.619 

2.743 

6.699 

1.353  -  0.823 

3.664 

2.462 

9.372 
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The  iron  region  has  cells  on  the  order  of  a  few  mean  free  paths  thick,  with 
the  exterior  water  tetrahedra  even  thicker.  The  errors  shown  in  Figure  88 
were  for  the  TETRAN  linear  characteristic  spatial  quadrature  for  each 
region.  In  Figure  89  to  Figure  91,  three  different  methods  are  compared  in 
each  region. 

All  the  methods  were  run  on  TETRAN  and  are  compared  to  the  MCNP 
solution  with  a  relative  uncertainty  as  shown.  A  spherical  harmonic  (SH) 
method  with  an  S8  quadrature  using  a  P3  expansion  was  run  using  the  LC 
method.  Note  that  SH  methods  can  create  negative  fluxes  independent  of  the 
spatial  quadrature  used.  When  this  happens,  the  EC  method  fails  to  run. 

The  other  two  methods  used  the  same  T-Scat  hbrary  employing  the  LC  and 
then  the  EC  method  for  direct  comparison  of  the  multigroup  anisotropic 
behavior. 


Figure  89.  Relative  difference  in  source  region  of  H20-Fe-H20  problem  to 

MCNP  solution. 
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Figure  90.  Relative  difference  in  middle  region  of  H20-Fe-H20  problem  to 

MCNP  solution. 


Outer  {H2O)  Region  Differences 


Figure  91.  Relative  difference  in  outer  region  of  H20-Fe-H20  problem  to 

MCNP  solution. 
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This  is  the  first  time  that  the  exponential  characteristic  method  has  been 
used  with  multi-group  anisotropic  scattering.  In  the  source  region,  the 
T-Scat  libraries  perform  well  in  the  first  group,  and  then  all  the  methods  are 
comparable.  In  the  iron  region  the  SH  method  is  markedly  better  in  the 
second  group  than  T-Scat  using  LC  or  EC.  Finally,  in  the  region  with  the 
thickest  tetrahedra,  both  LC  methods  are  comparable  while  the  EC  method  is 
about  ten  times  better  in  the  source  group  and  half  the  error  in  the  second 
group. 

It  is  expected  that  the  EC  method  would  perform  best  in  those  regions 
where  the  cells  are  thickest.  And  again  the  EC  method  can’t  be  used  with 
spherical  harmonics  because  of  the  creation  of  negative  sources.  The  result 
for  the  flux  in  group  1  (which  is  the  group  in  which  the  source  emits)  is  most 
like  a  deep  penetration  problem,  because  the  source  is  in  the  inner  box. 

Lower  groups  have  downscatter  sources  distributed  spatially  throughout  the 
problem.  Looking  at  Figure  91,  we  see  that  indeed  the  EC  method  performs 
much  better  than  the  other  methods  in  the  source  group.  Still,  one  might  be 
tempted  to  use  the  LC  method  and  tolerate  the  differences  in  error.  While 
the  EC  method  will  fail  if  a  negative  source  is  encountered,  the  LC  method 
will  include  it  in  its  and  continue  on.  The  end  result  can  be  a  negative  flux 
calculation,  which  is  nonsense.  This  is  demonstrated  with  the  next  problem. 

Again  the  triple  cube  geometry  is  employed.  This  time  the  materials  are 
water,  boron  (^°B)  and  water,  with  the  boron  density  being  10  g/cm^.  This 
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density  is  used  to  enhance  the  cell  optical  thickness  for  the  same  mesh  as 
used  earlier.  The  average  cell  optical  thicknesses  are  again  shown  in  Table  6. 


Table  6.  Average  cell  optical  thickness  for  H20-i®B-H20  system. 


Group  (MeV) 

Water/Source 

lOB  (10g/cm3) 

Water 

17.0  -  15.0 

.828 

9.852 

2.118 

15.0-13.5 

.882 

9.910 

2.255 

13.5-12.0 

10.00 

2.310 

12.0  -  10.0 

.946 

2.421 

10.0  -  7.79 

.970 

10.07 

2.480 

Compared  to  the  thicknesses  of  the  previous  problem,  the  water 
tetrahedra  are  thinner.  This  is  because  of  the  higher  energy  groups  used. 
This  increase  in  energy  also  increases  the  anisotropy  of  the  problem.  Figure 
92  shows  the  scalar  flux  for  the  problem  regions  found  with  T-Scat  LC  and 
EC  methods. 


For  the  water-iron-water  problem  the  scalar  flux  was  reduced  by  about 
three  orders  of  magnitude  through  the  shield  region.  Here,  the  reduction  is 
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five  orders  of  magnitude.  For  this  problem,  the  LC  method  produced 
negative  scalar  fluxes  in  the  outer  water  region  using  both  the  SH  and  T-Scat 
libraries.  Hence,  only  the  EC  method  provided  positive  results  for  the 
external  region.  And,  of  course,  EC  requires  a  strictly  non-negative 
scattering  matrix  provided  by  the  MC  facet  method  (SH  fails  in  this  regard). 

As  for  the  water-iron-water  problem  the  errors  are  shown  in  the  figures 
below.  Here,  the  LC  and  EC  methods  are  compared  except  in  the  outside 
water  region  where  the  LC  methods  failed  to  provide  positive  results. 

Instead,  a  finer  spatial  mesh  was  used  (1068  tetrahedra).  Again  the  LC 
methods  produced  negative  results,  but  the  overall  errors  in  the  outer  region 
were  reduced. 


Source  Region 


Figure  93.  Scalar  flux  differences  of  source  region  relative  to  MCNP 

for  the  H2O-10B-H2O  problem. 
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Figure  94.  Scalar  flux  differences  of  middle  region  relative  to  MCNP 

for  the  H20-1®B-H20  problem. 


Outer  Region 
H2O-B-H2O 


Figure  95.  Scalar  flux  differences  of  outer  region  relative  to  MCNP 

for  the  H20-1®B-H20  problem. 

While  the  errors  in  the  external  region  are  high  in  the  first  energy  groups, 
examination  of  Figure  92  shows  that  the  EC  method  has  done  an  excellent 
job  in  tracking  with  the  MCNP  solution  for  the  deep  penetration  problem 
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considering  the  substantial  attenuation  of  the  flux.  A  further  advantage  of 
the  discrete  ordinates  method  is  that  it  can  provide  a  profile  of  the  flux 
(dependent  on  the  refinement  of  the  mesh).  It  is  possible  to  take  a  closer  look 
at  exactly  where  the  flux  goes  negative  for  the  LC  problem  by  using  this 
ability. 

I  developed  a  simple  plotting  routine  that,  for  a  given  point  in  space, 
determines  the  closest  tetrahedron  center  and  reports  that  tetrahedron’s 
scalar  flux.  Figure  96  plots  the  group  1  scalar  flux  along  the  x-axis  passing 
through  the  center  of  the  problem.  The  charts  are  expanded  to  show  the 
negativity  just  as  the  exterior  water  region  is  entered.  The  EC  method 
maintains  positivity  for  this  problem  throughout. 


Figure  96.  Scalar  flux  map  of  the  tetrahedra  along  the  x-axis  of  the  water- 
boron-water  problem.  LC  method  fails  to  maintain  positivity. 


Plotting  the  EC  method  alone  on  a  log  scale  shows  the  moderating  effect  of 
the  water  in  the  external  region.  This  is  shown  in  Figure  97  for  the  first 
group  (15.0-17.0  MeV)  and  the  10‘h  group  (2.232-1.738  MeV). 
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Figure  97.  Scalar  flux  from  TETRAN’s  EC  method  for  the  water-boron-water 

problem. 

Returning  to  the  water-iron-water  problem,  we  apply  the  same  type  of 
chart  along  the  outside  surface  of  the  cube  to  get  a  profile  of  the  flux  at  the 
surface  (comparing  the  spherical  harmonic  and  T-Scat  methods).  Recall  that 
here  the  scalar  flux  was  positive  for  the  exterior  region.  Figure  98  shows  that 
at  the  corners  of  the  cube  the  scalar  flux  was  actually  negative.  Referring  to 
Figure  87,  the  flux  using  the  LC  method  was  reported  as  positive,  but  lower 
than  the  MCNP  result  in  the  external  region  across  all  groups. 
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Figure  98.  Scalar  flux  goes  negative  for  water-iron-water  problem  too. 


The  scalar  fluxes  reported  by  the  EC  method  are  not  reduced  by  the  non¬ 
physical  negative  flux  results  obtained  using  the  LC  method. 


Transport  Summary 

This  chapter  has  presented  only  a  representative  number  of  materials 
tested  with  scattering  matrices  created  using  the  MC  facet  method.  Its 
purpose  was  not  to  validate  the  TETRAN  code,  but  rather  to  demonstrate  the 
viability  of  the  non-negative  cross-section  libraries  we’ve  created.  We’ve 
shown  in  particular: 

1.  a  general  reduction  in  scalar  flux  errors  (as  compared  to  SH  methods)  for 
comparable  quadratures, 
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2.  particular  improvement  over  SH  methods  with  small  numbers  of 
ordinates, 

3.  ability  of  the  MCSN  code  to  use  T-Scat  multi-group  anisotropic  scattering 
data, 

4.  comparable  to  better  performance  (than  SH  methods)  for  fine  energy 
groups  (noting  that  ordinate-to-ordinate  methods  lack  angular  support  for 
such  group  structures), 

5.  the  ability  to  use  single,  mixed  and  multiple  materials,  and 

6.  the  ability  to  use  the  EC  method  on  deep  penetration  problems 

•  demonstrating  overall  superior  performance  to  LC/SH  methods  and 

•  guaranteeing  flux  positivity  with  non-negative  cross-section  matrices. 

These  observations,  together  with  those  of  Chapter  9  demonstrate  that 
the  MC  facet  method  is  a  viable  technique  for  general  use  in  discrete 
ordinates  (and  MCSN)  transport  codes. 
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11.  Conclusions 


Since  the  1970’s  researchers  in  the  transport  community  have  recognized 
the  problem  with  the  use  of  spherical  harmonics  to  approximate  the 
scattering  source.  Such  harmonics  can  create  negative  sources  which,  in 
turn,  create  non-physical  negative  fluxes.  This  research  has  shown  how  such 
fluxes  negatively  impact  the  calculation  of  the  scalar  flux.  In  particular, 
negative  fluxes  can  cause  catastrophic  failure  of  transport  codes  employing 
the  exponential  characteristic  (EC)  method.  Such  characteristic  methods 
have  only  been  developed  within  the  last  few  years  (here  and  at  the  national 
labs)  to  make  the  study  of  deep  penetration  problems  feasible  with  the 
discrete  ordinates  method.  The  inability  of  such  methods  to  perform  multi¬ 
group,  anisotropic  scattering  because  of  the  negativity  of  spherical  harmonics 
has  only  added  to  the  requirement  that  a  new  method  be  developed. 

A  number  of  attempts  have  been  made,  some  quite  recently,  to  solve  the 
problem  of  negative  scattering  cross-sections.  Unfortunately,  as  we  have 
seen,  each  has  its  major  drawbacks.  This  research  has  attempted  to  stay  true 
to  the  physics,  and,  in  so  doing,  has  avoided  the  need  for  complicated 
interpolations,  or  negative  fix-ups  to  correct  for  computational  complexities. 
The  result  is  a  method  that  has  significant  advantages  over  its  predecessors, 
is  simple  to  implement  and  provides  quality  results. 
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By  just  examining  the  scattering  matrices  themselves,  weVe  shown  that 
the  MC  facet  method  creates  a  T-Matrix  that 

•  is  consistent  with  the  physics, 

•  has  an  easily  adjustable  accuracy, 

•  guarantees  non-negative  cross  sections, 

•  guarantees  within  group  and  next  group  scatter,  and 

•  can  be  rebalanced  to  conserve  the  first  angular  moment  of  the  scatter. 
By  examining  the  use  of  the  T-Matrix  in  a  discrete  ordinates  transport  code 
we’ve  shown  that  the  method  can 

•  reduce  scalar  flux  errors  (as  compared  to  SH  methods)  particularly  for 
comparable  low  numbers  of  ordinates, 

•  facilitate  the  use  of  MCSN-like  codes  to  use  multi-group  anisotropic 
scattering  data, 

•  improve  performance  for  fine  energy  groups, 

•  easily  accommodate  single,  mixed  and  multiple  materials,  and 

•  be  used  by  transport  codes  employing  the  EC  method  for  deep 
penetration  problems. 

Overall  the  method  promises  to  be  viable  for  neutral  particle  transport  of 
all  types,  including  photon  transport.  Any  transport  scatter  that  can  be 
modeled  with  MCNP  can  be  modeled  using  the  MC  facet  method. 
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Appendix  A:  Sample  Mean  Monte  Carlo  Integration 

This  appendix  provides  a  brief  overview  of  sample  mean  MC 
integration.  This  technique  is  applied  to  the  average  scattering  integral  to 
calculate  the  T-matrix.  The  description  here  is  completely  general  and 
follows  the  description  of  Kalos^"^. 

We  wish  to  evaluate  an  integral  of  the  form 

G  =  {^g(x)fix)dx  (199) 

Ja 

where 

f(x)>0,  ff(x)dx  =  l.  (200) 

Ja 

Sample  the  probability  distribution  function,  f(x),  to  obtain  a  set  of 
variables  Xi,  X2,  X3...,  Xn,  and  form  the  arithmetic  mean 

=  G  +  (201) 

It  is  assumed  that  cr^  ,  and  Gm  are  normally  distributed  (by  the  law  of 

LjN 

large  numbers).  Thus,  the  uncertainty  estimate  can  be  specified  using 
confidence  intervals.  For  example,  there  is  a  90%  chance  that 

^G^|<1.65-^  (202) 

^JN 

where 

cr^  =  g^(x)f(x)dx-G^  .  (203) 

Ja 
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Reapplying  MC  integration  to  this  integral,  the  error  estimate  of  the  solution 
is 


<7jV  = 


— 

NT^ 


\2 

y 


il/2 


or 


(204) 


|^Gjv|<1.65 


—  Vrtf  _  —Yq. 

-sr^Sl 
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N‘ 


1/2 


4n 


(205) 


The  corresponding  multipher  for  a  95%  confidence  is  1.96,  and  for  a  99% 
confidence,  2.58. 

If  /(^does  not  appear  explicitly  in  the  integral,  one  can  simply  use 
f{x)  =  1/(6 -a)  and =  (6-a)  x  integrand.  This  choice  of  f{x)  requires  a 


uniform  random  sampling  in  [a,  6].  It  is  easily  implemented,  but  other 
choices  can  improve  efficiency  by  reducing  cr^ .  This  is  a  fundamental 
variance  reduction  scheme.  For  more  information  concerning  variance 
reduction  the  reader  is  referred  to  Kalos^®  or  Lewis  &  Miller^. 
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Appendix  B:  The  Inelastic  Neutron  Scattering  Laws  of  MCNP 


The  full  list  of  neutron  scattering  laws  employed  by  MCNP  follows. 

1.  Equi-probable  energy  bins  (Law  1) 

2.  Discrete  photon  energy  (Law  2) 

3.  Tabular  Distribution  (Laws  4  &  44) 

4.  General  evaporation  spectrum  (Law  5) 

5.  Simple  Maxwell  fission  spectrum  (Law  7) 

6.  Evaporation  Spectrum  (Law  9) 

7.  Energy  Dependent  Watt  Spectrum  (Law  11) 

8.  Tabular  linear  functions  (Law  22) 

9.  Equiprobable  energy  multipliers  (Law  24) 

10.  N-body  phase  space  distribution  (Law  66) 

11.  Correlated  energy-angle  scattering  (Law  67) 

Other  special  handling  techniques  include 

Fission  Inelastic  Scattering 
S(a,P)  Thermal  Scattering 

MCNP  also  includes  photon  and  charged  particle  interactions.  The  MC  facet 
approach  can  readily  be  applied  to  these  other  interactions. 
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Appendix  C:  Scattering  Levels  of  Oxygen  (^®0) 
Oxygen- 16:  Level  Inelastic  Energy  Loss  in  Collision 


Level 

Q 

n,  n1 

-6.0494 

n,  n2 

-6.12989 

n,  n3 

-6.9171 

n,  n4 

-7.11685 

n,  n5 

-8.8719 

n,  n6 

-9.585 

n,  n7 

-9.8445 

n,  n8 

-10.356 

n,  n9 

-10.957 

n,  n10 

-11.08 

n,  n11 

-11.0967 

n,  n12 

-11.26 

n,  n13 

-11.52 

n,  n14 

-11.6 

n,  n15 

-12.049 

n,  n16 

-12.44 

n,  n17 

-12.53 

n,  n18 

-12.796 

n,  n19 

-12.9686 

n,  n20 

-13.15 

n,  n21 

-13.45 

n,  n22 

-13.75 

n,  n23 

-14.05 

n,  n24 

-14.35 

n,  n25 

-14.65 

n,  n26 

-14.95 

n,  n27 

-15.25 

n,  n28 

-15.55 

n,  n29 

-15.85 

n,  n30 

-16.15 

n,  n31 

-16.45 

n,  n32 

-16.75 

n,  n33 

-17.05 

n,  n34 

-17.35 

Hi  n35 

-17.65 

n,  n36 

-17.95 

n,  n37 

-18.25 

n,  n38 

-18.55 
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Appendix  D:  LANL  30  Group  Structure 


Energies  for  the  LANL  30  group  structure  are  in  eV. 


1 .39E-04 

1.52E-01 

4.14E-01 

1.13E+00 

3.06E+00 

8.32E+00 

2.26E+01 

6.14E+01 

1.67E+02 

4.54E+02 

1.24E+03 

3.35E+03 

9.12E+03 

2.48E+04 

6.76E+04 

1.84E+05 

3.03E+05 

5.00E+05 

8.23E+05 

1.35E+06 

1.74E+06 

2.23E+06 

2.87E+06 

3.68E+06 

6.07E+06 

7.79E+06 

1  .OOE+07 

1.20E+07 

1.35E+07 

1.50E+07 
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Appendix  E:  VITAMIN-J  175  Group  Structure 


Energies  for  the  VITAMIN-J  group  structure  are  in  eV. 


1.00E-05 

1.00E-01 

1.13E+00 

3.93E+00 

1.37E+01 

4.79E+01 

1.67E+02 

5.83E+02 

2.03E+03 

3.04E+03 

7.10E+03 

1.93E+04 

2.61  E+04 
4.09E+04 
7.20E+04 
1.11E+05 
1 .43E+05 
1 .83E+05 
2.35E+05 
2.97E+05 
3.88E+05 
5.50E+05 
7.07E+05 
9.07E+05 
1 .22E+06 
1 .57E+06 
2.02E+06 
2.37E+06 
2.87E+06 
4.07E+06 
5.49E+06 
6.70E+06 

8.61  E+06 
1.11E+07 
1 .35E+07 
1 .57E+07 


4.14E-01 
1.45E+00 
5.04E+00 
1.76E+01 
6.14E+01 
2.14E+02 
7.49E+02 
2.25E+03 
3.35E+03 
9.12E+03 
2.19E+04 
2.70E+04 
4.63E+04 
7.95E+04 
1.17E+05 
1 .50E+05 
1 .93E+05 
2.47E+05 
2.98E+05 
4.08E+05 
5.78E+05 
7.43E+05 
9.62E+05 
1.29E+06 
1.65E+06 
2.12E+06 
2.39E+06 
3.01  E+06 
4.49E+06 
5.77E+06 
7.05E+06 
9.05E+06 
1.16E+07 
1.38E+07 
1.65E+07 


5.32E-01 
1.86E+00 
6.48E+00 
2.26E+01 
7.89E+01 
2.75E+02 

9.61  E+02 
2.49E+03 
3.71  E+03 
1 .03E+04 
2.36E+04 
2.85E+04 
5.25E+04 
8.25E+04 
1.23E+05 
1 .58E+05 
2.02E+05 
2.73E+05 
3.02E+05 
4.50E+05 
6.08E+05 
7.81  E+05 
1.00E+06 
1 .35E+06 
1.74E+06 
2.23E+06 
2.47E+06 
3.17E+06 
4.72E+06 
6.07E+06 
7.41  E+06 
9.51  E+06 
1.22E+07 
1 .42E+07 
1.69E+07 


6.83E-01 
2.38E+00 
8.32E+00 
2.90E+01 
1.01  E+02 
3.54E+02 
1 .23E+03 

2.61  E+03 

4.31  E+03 
1.17E+04 
2.42E+04 
3.18E+04 
5.66E+04 
8.65E+04 
1 .29E+05 
1 .66E+05 
2.13E+05 
2.87E+05 
3.34E+05 
4.98E+05 
6.39E+05 
8.21  E+05 
1.11  E+06 
1 .42E+06 
1 .83E+06 

2.31  E+06 
2.59E+06 
3.33E+06 
4.97E+06 
6.38E+06 
7.79E+06 
1  .OOE+07 
1.25E+07 
1.46E+07 
1.73E+07 


8.76E-01 
3.06E+00 
1.07E+01 
3.73E+01 
1 .30E+02 
4.54E+02 
1.58E+03 
2.75E+03 
5.53E+03 
1.50E+04 
2.48E+04 
3.43E+04 
6.74E+04 
9.80E+04 
1.36E+05 
1.74E+05 
2.24E+05 
2.95E+05 
3.69E+05 
5.23E+05 
6.72E+05 
8.63E+05 
1.16E+06 
1.50E+06 
1.92E+06 
2.35E+06 
2.73E+06 
3.68E+06 
5.22  E+06 
6.59E+06 
8.19E+06 
1.05E+07 
1.28E+07 
1.49E+07 
1.96E+07 
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